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We study the effects of focaf inhomogeneities, i.e., siow sites of hopping rate q < 1, in a totaliy 
asymmetric simpie exciusion process (TASEP) for particles of size £ > 1 (in units of the lattice 
spacing). We compare the simulation results of £ — 1 and £ > 1 and notice that the existence 
of local defects has qualitatively similar effects on the steady state. We focus on the stationary 
current as well as the density profiles. If there is only a single slow site in the system, we observe a 
significant dependence of the current on the location of the slow site for both 1=1 and I > 1 cases. 
When two slow sites are introduced, more intriguing phenomena emerge, e.g., dramatic decreases in 
the current when the two are close together. In addition, we study the asymptotic behavior when 
q —* 0. We also explore the associated density profiles and compare our findings to an earlier study 
using a simple mean- field theory. We then outline the biological significance of these effects. 
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I. INTRODUCTION 



A better understanding of non-equilibrium steady 
states in interacting complex systems forms a critical goal 
of much current research in statistical physics. In this 
pursuit, the totally asymmetric simple exclusion process 
(TASEP) [T], [2, y , H, \M, m nas played a paradigmatic role. 
It provides a nontrivial, yet exactly solvable, example of 
phase transitions far from equilibrium, taking place even 
in one-dimensional (ID) lattices. At the same time, it 
also serves as the starting point for the modeling of many 
physical (driven diffusive ^processes, such as translation 
in protein synthesis inhomogeneous growth pro- 

cesses (e.g. Kardar-Parisi-Zhang growth) [ToLTlll] and ve- 
hicular traffic 0, [H . 

In its simplest version, the TASEP involves a single 
species of particles hopping to nearest-neighbor sites, in 
one direction only, along a homogeneous ID lattice. Pro- 
vided the destination site is empty, the rate for the par- 
ticle hop is fixed at 7 (typically chosen as unity with- 
out loss of generality). With periodic boundary condi- 
tions, the steady-state distribution is trivial [l4[ but the 
full dynamics is quite complex [lj| [H, G3 • With open 
boundary conditions, particles are injected with rate a 
(in units of 7) at one end and drained with rate [3 at 
the other end. The competition of injection, transport 
and drainage induces a nontrivial phase diagram in the 
a- (3 plane [l], [3, H, 0, H, @| , reflecting a highly nontrivial 
steady state. Three phases are present: a maximum- 
current phase for a, (3 > 1/2, and a low- (high-) density 
phase for a < 0, a < 1/2 ((3 < a, (3 < 1/2). Not surpris- 
ingly, there are also rich dynamical aspects [H, [n| . 

To model protein synthesis, each site on the lattice 
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represents a codon on the messenger RNA (mRNA), and 
the particles represent the ribosomes. Injection, hop- 
ping, and drainage are associated respectively with ini- 
tiation, elongation, and termination in biological terms. 
The quantity of interest, namely, the (steady-state) pro- 
tein production rate, is identical to the (stationary) par- 
ticle current. Clearly, the simple TASEP falls short of 
the biological system in several significant aspects. One 
is that an individual ribosome "covers" several codons 
@j HH) as opposed to a particle occupying only a 
single site. Another is that, in all naturally occurring 
mRNAs, the codons carry genetic information and there- 
fore necessarily form an inhomogeneous sequence. Thus, 
the elongation rate of a ribosome is unlikely to be uni- 
form; instead, the hopping rate, 7$, of a particle becomes 
a function of the site i. For example, it is well known 
that translation slows down at specific codons (see, e.g. 
d S3, [H, H HI), with potentially significant conse- 
quences for protein production rates. Indeed, the steady- 
state current may depend sensitively on not only the fre- 
quency of each codon's occurrence, but also the order of 
their appearance in the sequence. Both of these issues - 
extended objects and inhomogeneous rates - have been 
addressed recently in separate contexts which we sum- 
marize briefly in the following. 

The results associated with inhomogeneous (quenched 
random) rates fall into two broad categories, in the sense 
that the randomness can be associated with the particles 
[H, H3i [HI or with the sites. Randomness of the for- 
mer type is more relevant for vehicular traffic where it 
accounts for a variety of driver preferences. In contrast, 
the disorder in the protein case is clearly site- dependent, 
leading to spatially non- uniform hopping rates 7^. Re- 
stricting ourselves to this class, we can consider the ef- 
fect of having a whole distribution, or very specific con- 
figurations, of 7i. Starting from given distributions, two 
groups [29|, H3| studied the resulting disorder-average in 
periodic systems. To mention just one significant effect, 
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the current-density diagram develops a plateau: limited 
by the smallest rate in the system, the current becomes 
independent of density over a range of densities. Harris 
and Stinchcombe [3(| also extended this work to open 
systems. While these studies may be of some interest to 
mixtures of many different mRNAs, our primary interest 
here is to understand how the production rate of a spe- 
cific protein is associated with a specific genetic sequence. 
As a first step towards a solution, we adopt the approach 
of several other studies i, [H H HI , by focusing on the 
effects of a few localized inhomogencitics, i.e., hopping 
rates which are uniform except at a handful of sites [3J] . 
As a synthesis of these studies, we will explore in some 
detail the consequences of having extended objects and 
locating one or two slow sites at a variety of positions 
on the lattice. In this manner, by introducing more and 
more sites with a range of rates, we hope to understand 
inhomogeneities in a systematic way, setting the stage for 
further investigations of the translation process. 

A full comprehension of the effects of slow sites on 
the particle current may have potentially significant ap- 
plications in biotechnology. While there are 64 distinct 
codons, proteins are chains composed of just 20 amino 
acids. So, many different mRNAs (codon sequences) 
can code for the same particular protein. Moreover, the 
amino acid is incorporated into the growing chain by 
an important intermediary, the so-called transfer RNA 
(tRNA), which carries the complementary anticodon. It 
turns out that the mapping between codons and tR- 
NAs is also not precisely 1-1. For example, in E. coli, 
the genetic code actually involves 61 sense codons and 
about 46 tRNAs with associated anticodons (3f| ■ Mean- 
while, for a given mRNA sequence, the protein produc- 
tion rate is often modeled in terms of (generally accepted) 
chargcd-tRNA (aminoacyl-tRNA, or aa-tRNA) concen- 
trations [HI, so that different sequences can result in 
different production rates for the same protein. By elu- 
cidating how the spatial distribution of defects, especially 
of bottlenecks, affects translation rates, we can pinpoint 
those clusters of codons which are likely to have the most 
significant effect on the production rate of the associated 
protein. Exploiting the degeneracy in the mapping from 
mRNA sequence to protein, we can provide guidance as 
to how a few selected, local modifications of the mRNA 
can optimize the production rate of a given protein. 

Our paper is organized as follows. In Section [Til we 
define the model and provide a more detailed description 
of previous work, concerning exclusion processes with ex- 
tended objects or spatially inhomogeneous rates. A brief 
discussion of a previous mean-field analysis is also in- 
cluded. In Section IIIH we present our Monte Carlo re- 
sults. We focus especially on the implications of hav- 
ing extended objects by varying the particle size. We 
first consider the interaction between one slow site and 
the system boundary, and, motivated by the resulting 
findings, turn to the interactions between two slow sites. 
This provides new insights for genes containing clusters 
of slow codons, which occur frequently in, e.g., E. coli, 



Drosophila, yeast and primates [9|, [36|, [37[ . In Section HVl 
a complete investigation of systems with inhomogeneities 
is presented using a mean-field approach. Section [V] con- 
tains our conclusions and a summary of open questions. 

II. MODEL SPECIFICATIONS AND KNOWN 
RESULTS 

The TASEP is defined on a ID lattice of N sites. We 
introduce an index i = 1,2, N to label the sites. Each 
site (codon) is either occupied by a single particle (ribo- 
some) of length I (in units of sites) or empty. A micro- 
scopic configuration of the system can be uniquely char- 
acterized in terms of a set of occupation variables, {n.i}, 
taking the value 1(0) if site i is occupied (empty). Of 
course, the extended nature of the particle induces strong 
correlations in {n,;}, in the sense that a single ribosome 
always covers t consecutive sites. Yet, at any given time, 
only one of the covered codons is being "read" (i.e., the 
codon is "covered by the aminoacyl site" , or A site, of 
the ribosome) and translated into an amino acid. Here, 
we refer to the associated location on the ribosome as 
the "reader" (of the genetic code). For our purposes, it 
is not essential which one of the I sites is labeled as the 
reader, and so we follow the convention in [8[ and choose 
the first (leftmost) site. Hence, the statement "a ribo- 
some (or particle) is located at site i" implies that the 
reader is located at site i and the subsequent 

I=£-l (1) 

sites are also occupied. Naturally, the position of the 
reader determines the elongation rate, i.e., ji, since the 
ribosome must wait for the arrival of the aa-tRNA with 
the i-specific anticodon before it can move to the next 
site. Clearly, the reader locations can also be used to 
label a microscopic configuration, i.e., we can define the 
reader occupation number at site i as r^. The sets {rii} 
and {ri} are uniquely related to each other. Moreover, 
due to the extended size of a particle, strict constraints 
are built in (e.g., = 1 implies 7'; + i = ... = rj + £_i = 
and rii = ... = ri;+£-i = 1). As a consequence, neither 
set can be arbitrary and serious correlations arise as soon 
as^>l[l|. 

In our simulations, we adopt a random sequential up- 
dating scheme and keep a list of locations of readers. In 
addition, the site i — is always occupied by a "virtual 
reader," which accounts for particles entering the system 
(initiation). At the beginning of each Monte Carlo step 
(MCS) , we first find the number of particles in the sys- 
tem and label it M . Then, we randomly select an entry 
from this list of M + 1 readers. If the chosen reader is 
virtual (i.e., i = 0), a new particle enters the lattice with 
probability a, provided all the first £ sites are empty. If 
the chosen reader is real, say, at site i > 0, the associ- 
ated particle is then moved to site i + 1 with probability 
7i, provided site i + £ is empty. With this notation, we 
can also write the initiation and termination probabilities 
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FIG. 1: Sketch of a TASEP for particle size £ = 6 with (a) a 
single slow site at position k, with rate q, and (b) two slow 
sites with rate q, separated by a distance d. 



(a and (3) as 70 and jn, respectively. To be complete, 
the sites beyond the lattice are by definition empty, so 
that once a particle reaches N — I + 1, it will not ex- 
perience steric hindrance (see Fig. Q] for a sketch of this 
process). These processes have been termed "complete 
entry" and "incremental exit" [39| . Other entry and exit 
rules can be considered, but are believed to be incon- 
sequential provided t/N <C 1. Each MCS consists of 
M + 1 such attempts, giving an even chance, on average, 
for each particle (ribosome) in the system to elongate or 
terminate, as well as for an initiation event to occur. 

Starting with an empty lattice, we typically discard 
2 x 10 6 MCS to ensure that the system has reached the 
steady state. Unless otherwise noted, good statistics re- 
sult if we average over least 2 x 10 4 measurements, sep- 
arated by 100 MCS in order to avoid temporal corre- 
lations. Such steady state averages will be denoted by 
(...). To reduce the number of parameters in the model, 
we study systems with a = ft = 7j = 1, except at one or 
two sites. The system sizes (N) range from 200 to 1000, 
with most data taken from N = 1000. 

To characterize the state of the system, we monitor 
several observables. The most obvious is p\ = (r*j), a 
quantity we will refer to as the ribosome (or reader, or 
particle) density. Of course, Yli Pi ^ s J us t the average 
number of particles in the system (i.e., ribosomes on the 
mRNA). Thus, the overall particle density N~ 1 J2iPi 
is bounded above by l/£. Another interesting variable 
Pi = (rii), labeled as the "coverage density", is the prob- 
ability that site i is covered by a particle (regardless of 
the location of the reader). Needless to say, the pro- 
file for the vacancies is given by the local hole density, 
p' 1 = 1 — pi. The overall coverage density, N^Y^iPu 
may reach unity and provides a good indication of how 
packed the system is. The two profiles are related by 
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A quantity of great importance to a biological system 
is the steady-state level of a given protein. If we assume 
that the degradation rates are (approximately) constant, 
i.e., independent of protein concentration, then these lev- 
els are directly related to the protein production rates. In 
our model, such a rate is just the average particle current 
J, defined as the average number of particles exiting the 
system per unit time. In the steady state, it is also the 
current measured across any section of the lattice. For 
simplicity and to ensure the best statistics, we count the 
total number of particles which enter the lattice over the 
entire measurement period (at least 2 x 10 6 MCS in most 
cases). 

In this study, we focus on two simple types of inhomo- 
geneities: one or two "slow" sites (Fig.[IJ. Their locations 
specify the only inhomogeneities in the rates. 

One slow site, at position k. We denote 7^ by q (< 1). 
This corresponds to a bottleneck in the lattice. We are 
especially interested in the dependence of the current, 
denoted by J q (k), on the parameters q and k. 

Two slow sites, at positions hi and hi with separation 
d = (&2 — fei) and rates qi.2 = 7^ 2 . We find that, when 
qi 7^ 92 , the current is controlled mainly by the smaller of 
the two, with little dependence on d, in agreement with 
a simple mean- field theory to be discussed in Section [TV] 
Therefore, most of our attention will be devoted to the 
case with q± = qi = q < 1. Moreover, we choose to limit 
our study to both sites being far from the boundaries. 
Then, the current is insensitive to their average position 
(fc2+fci)/2, and we can investigate J q {d). Note that these 
are precisely the systems studied in Q, except that we 
consider particles with a range of sizes: £=1,2, 4, 6, and 
12. While there are qualitative similarities, we will dis- 
cuss the quantitative differences due to £ > 1, as well as 
interesting phenomena associated with the density pro- 
files. 

Let us now provide the context of our work by briefly 
reviewing some related earlier studies. The homoge- 
neous case (71 = ... = 7jv_i = 1) with I = 1 is ex- 
actly soluble @, H, 0, H, @] , and displays three phases in 
the a- (3 phase diagram. For I > 1, no exact solutions 
exist [40(. Analytic approximations using various mean- 
field approaches @, H HO, [H, Hl[ predict the presence of 
the same phases, though the phase boundaries depend 
on I (Fig. [2]) through the combination Q 
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(3) 



with the understanding p\ = for i < 0. 



Monte Carlo studies [8|, l39( largely confirm these conclu- 
sions. 

The three phases carry different currents and display 
distinct density profiles [1, S 0, S 0, 1, [H [H . Apart 
from "tails" near the boundaries, the (coverage) den- 
sity profiles approach uniform bulk values in the ther- 
modynamic limit, i.e., pi — > Pbuik, for 1 <C i <C N. 
For a < x an d a < (3, the system is in a low-density 
phase (L), characterized by pbuik = &*/ \X + ax) and 
J = a(l — a)/ (l + at). A high-density phase (H) pre- 
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FIG. 2: Phase diagram for an ordinary TASEP. On the dashed 
line, the H and L phases coexist. 



TABLE I: J-ptuik relation for particles of size I (£ : 



1). 



phase 



current J 



bulk density pbuik 



L 
H 

M 



q(1 - a) I (1 + at) 
13(1- (3)/ (1 + /31) 



■x 2 



la /(l + al) 
1-/3 

i-x 



vails for j3 < x an d j3 < a, with bulk density pbuik = 1 — 
and current J = [3(1 - (3)/ {l + (3t). For a, (3 > x, 
the system is in a maximum-current phase (M), where 
Pbuik = I — X and J — x 2 - On the a = (3 < x lme 
(dashed line in Fig. [5]), the system consists of two macro- 
scopic regions, characterized by a low (high) density re- 
gion near the entry (exit) point. The two regions are 
joined by a shock front that performs a random walk. 
This is often referred to as the "shock phase" (S). Ta- 
ble U summarizes the J — pbuik relation for TASEP with 
extended objects. There is good agreement between sim- 
ulations (with £ < 12) and analytic results for these bulk 
quantities [1, The details of the profile for i > 1, 
especially near the lattice boundaries, are less well un- 
derstood. While periodic structures (of period I) can be 
expected, mean-field theories @, H(| were successful 
in capturing only a limited part of the phenomena ob- 
served. We will return to these considerations in Section 
IIVI Beyond homogeneous systems, several studies intro- 
duced one or more "impurities" into TASEP with peri- 
odic boundary condition. A single "slow" site induces a 
shock in the densit y p rofile with some interesting statis- 
tics [H, 0, 145L w\ . Subsequently, generalizations to 
systems with a finite fraction of slow sites, randomly lo- 



cated, were also investigated [29]. For the richer case of 
the open boundary TASEP % M, [S3, EH 51 > Kolomeisky 
focused on point particles (£ = 1), with a single impu- 
rity at the center of the lattice [3l|, so as to mimic a 
defect situated deep in an infinitely long system. The 
consequences of the defect having both faster (q > 1) 
and slower (q < 1) rates were explored. By matching 
two ordinary TASEPs across the defect, the properties 
of such systems in the a-(3 plane can be well described 
[3ll . l48j |. While a fast site has no effect on the phase 
diagram, a slow site leads to a shift of the M-H and M- 
L phase boundaries to ^-dependent, smaller values of a 
and [3. The density profiles are quite sensitive to the ex- 
istence of a defect site [H EH). Kolomeisky's approach 
was generalized to the I = 12 case in (33|, with similar 
levels of success. Below, we will provide further details 
of this work, on which we base much of the analysis of 
our problem. Ha and den Nijs also studied the £ — 1 
open boundary TASEP with a single defect at the center 
[32l |. Focusing on the multicritical point a = /? = 1/2, 
they were mainly interested in the so-called "qucueing 
transition" and its critical properties. Detailed results of 
density profiles, such as power law behavior and critical 
exponents, were obtained in the region q = q c . Here, 
q c denotes the critical value of q below which the bulk 
density in front of the slow site deviates from the den- 
sity behind the blockage. By contrast, our focus here is 
essentially that of 0, [48[ , namely, how does the number 
and the locations or spacings of the slow sites affect the 
current through the system? In the single slow site case, 
investigated mainly the overall current as a function 
of q, while [48] analyzed the current as a function of k 
(the distance between the defect and the entry point). 
For the case of two defects, both studies find that the 
spacing between them plays a significant role for the cur- 
rent. A finite-segment mean- field theory in Q provides 
excellent agreement with data. In particular, clustered 
defects reduce the current much more effectively than 
well separated ones. In this sense, we can regard these in- 
vestigations as exploring the "interactions" between the 
slow site(s) and/or the boundaries. Since both studies 
are restricted to point particles (£ = 1), our intent is 
to explore the effects of having extended objects (with 
£ < 12). Though we expect qualitatively similar behav- 
ior, as pointed out in we also find noteworthy quanti- 
tative differences. Finally, we should mention that larger 
numbers of defects do not lead to significantly different 
effects Q, so that we limit ourselves to one or two slow 
sites here. 



III. MONTE CARLO RESULTS 

In this section, we present our Monte Carlo results. 
For convenience, we use a consistent color coding scheme 
(online only) for the various particle sizes, as specified in 
Table |TTJ The data here consist of the overall currents 
and the density profiles of both coverage and ribosomes. 



5 



TABLE II: Color coding scheme 



size £ 


online color 


1 


black 


2 


red 


4 


brown 


6 


green 


12 


blue 



Our focus will be how these quantities depend on q, k (for 
the case with a single defect), and d (for the case with 
two defects) . Although the profiles arc difficult to extract 
experimentally, the reader profiles will be of interest in 
subsequent studies involving real gene sequences, since 
they provide information on how frequently a particular 
tRNA is bound to the mRNA. By contrast, the currents 
are easily measurable and our results here may generate 
more immediate interest. 



A. One slow site 

We begin by placing one slow site on the lattice as in 
Fig. [Ha). Fig. [3] shows several coverage density profiles 
for a typical choice of parameters: N = 1000, q = 0.2, 
and k = 82 with £ = 1,6, 12. As expected, we observe 
pile-ups of particles due to the blockage - a high (low) 
density region before (after) the bottleneck in all three 
cases. However, due to the lack of ordinary particle-hole 
symmetry in the £ > 1 cases, the average densities on 
either side of the slow site are not symmetric around 
0.5. Instead, they are roughly related through the J-pbuik 
relations in the H and L phases, summarized in Table [J 
In detail, the profiles are quite different: The "tails" , i.e., 
the deviations from the bulk values, are quite noticeable 
in the vicinity of both the slow site and the edges of the 
system for the £ — 6,12 cases. The inset exposes more 
clearly that there are period £ structures in the profiles 
0, HojH^], especially just before the slow site. 

A more dramatic difference between point particles and 
extended objects emerges when we plot the ribosome den- 
sity p\, in Fig. 21 corresponding to the inset in Fig. [3] 
Similar to profiles in @, H3, , we find distinct period £ 
structures before the slow site. While the reader "waits" 
to pass the blockage, the readers of the following parti- 
cles tend to catch up and pause at sites k — n£, where 
n = 1, 2, ... The "tails" are even more marked than those 
in Fig.[3J To emphasize the difference between the reader 
and coverage profiles (p\ and pi), we show a case with 
q = 0.05, k = 948, £ = 12 in Fig. [5J Though both profiles 
contain the same information, we see that p\ (lower plot) 
is far more sensitive than pi (upper plot) in showing the 
very long tails (~ 1000 in this example) hidden in the col- 
lective behavior of the particles. At present, the crucial 
ingredients that control the characteristic decay length of 
the pJJ-envelopes have not yet been identified. Certainly, 
these very large length scales are completely absent from 
the £ = 1 systems deep within the H/L phases. 
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FIG. 3: (Color online) Coverage density profiles with one slow 
site of q = 0.2 at k = 82. I = 1,6, 12 (from bottom to top 
in both subsections) and N = 1000. The inset is a magni- 
fied view of the i £ [1, 150] interval, to expose the period £ 
structures. 




FIG. 4: (Color online) Ribosome density profiles with one 
slow site of q = 0.2 at k = 82. t = 1,6,12 (from top to 
bottom in both subsections) and N = 1000. Only the first 
150 lattice sites are shown. 



As for the current, Fig. [6] illustrates its dependence on 
q, fc, and £. Not surprisingly, the current is limited by 
the bottleneck and therefore varies monotonically with q. 
It is also reduced if the particle size increases, an effect 
that can be traced mainly to the particle density being 
effectively lower by the factor £. For point particles, the 
current is not very sensitive to the location of the slow 
site. The enhancement as k approaches the boundary 
of the system - referred to as the "edge effect" [48[ - is 
quite small. For larger £, the enhancement is much more 
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FIG. 5: (Color online) Coverage density profile (top) and ri- 
bosome density profile (bottom) with one slow site of q — 0.05 
at k = 948. I = 12 and N = 1000. 
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FIG. 6: (Color online) J q (k) as a function of the location k 
of the slow site for q — 0.2 (lower set of squares); 0.3 (middle 
circles) and 0.4 (upper triangles), (a) I = 1; (b) I — 2; (c) 
I = 6; (d) I = 12. In all cases, JV = 1000. 



pronounced, especially for smaller q. Whatever the mag- 
nitude, in all cases the current increases monotonically 
as the slow site is located closer and closer to the entry 
point. For 1=1, particle- hole symmetry is manifest in 
the microscopic dynamics, so that the symmetry of J q (k) 
under k — ■> N + 1 — k inversion, is obvious [48| . For £ > 1, 
the density profiles confirm the lack of this particle-hole 
symmetry very clearly. Correspondingly, there is a sys- 
tematic asymmetry in the current: J q {k) = J q (N + 1 — k) 
is satisfied only for k < I. The origin of this behavior is 
not well understood. 

The edge effect, and specifically its dependence on q 



FIG. 7: (Color online) Ai(q) for I = 1,2,4,6 and 12. 



and £, can be quantified by the ratio: 



Ai(?) = 



J q (k 



1) 



J q (k — > OO) 



(4) 



Fig. [7] shows that A-i(q) depends on q in a nontrivial 
way. The maxima of Ai(g) occur at lower values of q 
as £ increases, reminiscent of the behavior of the phase 
boundary between M and L/H. With appropriate scaling, 
the curves of Ai(g) can be collapsed for large £'s. From 
the biological perspective, the edge effect is not easily ob- 
servable since the current enhancement is less than 10% 
for the relevant I. 

Returning to Fig. [SJ we note that significant deviations 
from the asymptotic value, J ? (oo), are found only for fc's 
within a small distance from the boundaries. At a ca- 
sual glance, this range appears to depend on both q and 
£. On closer examination of, say, the most prominent 
case here: (q,£) = (0.2,12), we find that the decay of 
J q {k) into J q {oo) fits an exponential quite well (Fig. [8]), 
i.e., J q (k) — Jg(oo) oc exp (— k/5), with 6 s» 10. Assuming 
this behavior persists in the other cases, we can study the 
(q, £) dependence of this characteristic length and denote 
it by S(q,£). We believe that the origin of this length 
scale can be traced to the presence of tails in the density 
profiles near the lattice boundaries: As the slow site ap- 
proaches the entry point, these tails may seriously affect 
the injection process, and thus, the current as well. For 
the £ = 1 case, we observe that 5(q,£) ~ £(/2e//)j where 
£ is the characteristic length associated with the bound- 
ary layer of the density profile in the ordinary TASEP. 
Specifically, with entry rate a = 1 and exit rate /?, the 
system is in the H-phase and the profile decays expo- 
nentially into the bulk, as p x — pbuik oc exp(— x/£) with 
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FIG. 8: (Color online) Dependence of J q (k) on k obtained 
from simulation is plotted in squares and the line is a linear 
fit with slope equals -0.11. q = 0.2, £ = 12 and N = 1000. 



Here, the left half of our system is such a TASEP, except 
that we have an effective (3: /3 e ff = q/(l + q) [33|. Us- 
ing these arguments on the three g's shown, we estimate 
decay lengths of about 5 (q = 0.4), 3 (q = 0.3), and 2 
(q = 0.2) lattice constants. Though the data on the dif- 
ferences Jq(k) — Jq(oo) are small and noisy, simulation 
results are consistent with S(q, 1) ~ £.(Peff)- However, 
for £ > 1, there is no analytic result for the boundary 
layers of the density profiles. Moreover, the data sug- 
gest that they are quite complex (e.g., in Figs. [4] and [5]). 
Thus, it is unclear how to quantify the picture for point 
particles to the general case of S(q, £). At present, a com- 
plete understanding of both "boundary layers" - in the 
density profiles and in J q (k) - remains elusive. 

If we consider the edge effect as an "interaction" be- 
tween the slow site and the lattice boundaries, the nat- 
ural next step is to explore the interactions between two 
slow sites. In order to avoid edge effects, we place the 
two slow sites sufficiently far away from the boundaries 
and vary their separation. 



B. Two slow sites 

As mentioned in the previous section, the currents in 
the qi q2 cases are essentially controlled by the slower 
of the two rates and so, may be regarded as systems with 
a single slow site. These profiles can be interesting, but 
we choose to restrict our attention here to a study of the 
qx = qi = q case, in which the currents show a nontrivial 
dependence on d, the distance between the two slow sites. 
With two bottlenecks, the system consists of three sec- 
tions: before the first blockage, in between the two, and 
after the second defect. Of course, for small q, the overall 
density in the first (last) section is expected to be high 
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FIG. 9: (Color online) Coverage density profiles for two slow 
sites with q = 0.2. I = 12, d = 100; £ = 6, d = 125; 1 = 2, 
d — 150; and I — 1, d = 170 (Curves are plotted from left to 
right in the mid-subsection and top to bottom elsewhere). In 
all cases, N = 1000. Inset: q = 0.6, I = 12, d = 2400 and 
TV = 7200. 



(low). In these cases, the effective entry and exit rates for 
the central section are also low, so that a wandering shock 
should be present. Hence, the average profile should be 
linear for £ = 1 (and essentially so for larger I ||) with 
a positive slope. This behavior is understandable, since 
the section between the two defects is comparable to an 
ordinary TASEP with small a — (3. These expectations 
are generally confirmed by simulations with q < 0.5 and 
various d's up to 12. Fig. [5] shows typical coverage pro- 
files, for a relatively small rate of q = 0.2. The system 
appears to make a transition from this H/S/L phase to 
an M/M/M phase as q increases. The center profiles be- 
come essentially flat, as illustrated in the inset (where 
q = 0.6 and £ — 12). Details of this transition are being 
explored. 

More interesting are the finer features of the profiles 
in the small q cases. As in the single defect system, the 
profiles exhibit period £ structures near the slow sites. 
To resolve these more clearly, we plot the reader density 
profiles in Fig. [TO] In all cases that involve extended par- 
ticles (£ > 1), the readers clearly pile up behind the slow 
sites. Apart from these "jams," another feature emerges, 
namely, a sequence of depletion zones, each of which pre- 
cedes one of the period £ peaks. For £ = 2, the differences 
between the upper and the lower envelope are especially 
dramatic. More remarkably, when the blockages are sep- 
arated by small d's, two different, "overlapping tails" are 
created, as illustrated in the inset of Fig. [101 where d = 1, 
q = 0.2, and £ = 12. Indeed, there are further interesting 
structures for d < £, which will be presented elsewhere. 

Compared to these remarkable characteristics in the 
profiles, the behavior of the currents seems lackluster. In 
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FIG. 10: (Color online) Ribosome density profiles with two 
slow sites of q = 0.2. I = 12, d = 100; t = 6, d = 125; t = 2, 
d = 150; and I = 1, d = 170 (Curves are plotted from left 
to right in the mid-subsection and bottom to top elsewhere). 
Inset, 1 = 2 and d = 1. In all cases, N = 1000. 



Fig. [TTJ we plot four sets of currents (5TJ, J q (d), associ- 
ated with £ = 1,2,6, and 12. In all cases, we see that 
J is considerably suppressed when d is reduced. When 
the slow sites are very far apart, the current behaves as if 
there is only one slow site, consistent with expectations 
from mean-field theories. At the other extreme, when 
the two defect sites are nearest neighbors, the current 
reaches its minimum. Not surprisingly, period £ struc- 
tures emerge as d is varied, illustrated in the inset of 
Fig- flTT d). but become less prominent for d > 50. These 
plots also reveal that, unlike the dependence on k above, 
there are serious deviations from the d — > oo values when 
d is decreased. To quantify this deviation, we define 



A 2 (g) 



J q (d = 1) 

J q (d — * OO) 



(G) 



and plot this quantity vs. q in Fig. 1121 In contrast to 
Ai(<7), we observe that A2(<7) exhibits a sizable depen- 
dence on q, especially for small values of q. In the limit 
of q — > the current decreases by a factor of 2! In the 
following section, we will see that this factor can be un- 
derstood via a mean-field approach. 

To summarize our simulation results, two bottlenecks 
near each other have a dramatic effect on the current. We 
may regard this phenomenon as an "interaction" between 
the two slow sites, inducing far more "resistance" when 
they are close than when they are well separated. 

Two additional comments are in order. First, we re- 
turn to one of the predictions of the mean-field theory, 
namely that a second slow site, spaced far apart from its 
partner, should have no further effect on the current. Our 
data indicate that the current for two slow sites, spaced 
far apart, is systematically lower than the current for a 
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FIG. 11: (Color online) J q (d) as a function of the separation d 
between the two slow sites for q = 0.2 (lower set of squares); 
0.3 (middle circles) and 0.4 (upper triangles), (a) I = 1; 
(b) I = 2; (c) I = 6; (d) £ = 12. The inset in (d) is a 
magnified view of the d £ [1, 60] interval, to expose the period 
£ structures. In all cases, N = 1000. 
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FIG. 12: (Color online) A 2 (<?) for £ = 1,2,4,6 and 12. 



single slow site, but only by a very small amount (less 
than 1%). Second, we can again attempt to identify a 
length scale which controls how J q {d) approaches J q (oo), 
as d increases. Since the central section of the system 
displays a shock, it is natural to ask whether the intrin- 
sic width of the shock sets this length scale. According to 
[42L [loT ] , this width covers only a few lattice spacings in 
the periodic TASEP with a single defect. Here, however, 
it appears that the shock is much broader. For example, 
the averaged profile of the shock for the case of q = 0.2 
with point particles is shown in Fig. 1131 as well as a sim- 
ple fit using a tanh function [46] with width of about 10. 
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FIG. 13: (Color online) The open circles mark the average 
profile of a shock between sites 149 and 349 with q = 0.2, 
compiled from of a very long run (3 x 10 s MCS) with an 
N — 1000, 1 = 1 system. Details of how raw profiles are 
shifted (so that the shock is located at site x = shown 
here) will be published elsewhere [49j]. A simple fit using 
A + B tanh (x/10) is also shown: solid line. 



Intriguingly, this length appears to be comparable to the 
one appearing in Fig. fTTT a). More work is needed to fully 
explore these issues. 



IV. MEAN-FIELD THEORETIC APPROACHES 

Mean-field theory is known to agree well with the ex- 
act results for a number of macroscopic quantities in 
the steady state of the I = 1 TASEP (see, for example 
Q). For extended particles, no exact solution is avail- 
able [13] so that mean-field (and more sophisticated clus- 
ter) approximations form the only route toward some un- 
derstanding of the system's behavior. However, there are 
many levels of "mean-field" approximations [ll fl US Ull j 
corresponding to neglecting different types of correla- 
tions. For certain quantities (e.g., currents in large sys- 
tems), predictions from the simplest level are very close 
to the simulation results. For others (e.g., some reader 
profiles), only the most sophisticated level performs ad- 
equately. In all cases, no level of mean-field theory can 
give a good fit to both the current and the profile. A 
thorough discussion of all of these schemes is quite in- 
volved and will be provided elsewhere (49l |. Here, we will 
restrict ourselves to the simplest method and compare its 
performance to the simulations. 

All approaches start with the exact expressions for the 
current 

J = a(l-m) (7) 
= H (n (1 - n i+i )) ; i€[l,N-£] (8) 



it [n) 



i € 



[N-e,N-i] 



(9) 
(10) 



In the absence of the steady-state distribution, the most 
naive approximation is to replace (rifij) by (r^) (rij). Un- 
fortunately, the constraints due to particles with I > 1 
are so severe that this approximation is entirely inad- 
equate when j = i + i. Even for the simple case of 
TASEP on a periodic ring, it leads to an erroneous ex- 
pression for J (except if I = 1). Instead, the average 
(coverage) density at site % + I is much larger than the 
(conditional) probability that it is actually covered given 
that the reader is at site i. MacDonald and Gibbs (MG) 
proposed 0] a much better approximation: 



(n (1 - n l+e )} 



Pii 1 - Pi+e) 



P\p\+i 



1 



p l+ i + p\ +l p\ +l + p\ +t 



As discussed in Section HH the densities far from bound- 
aries are uniform (e.g., pJL^oo — > Pbulk /£) and this fact 
provides a good description of the current-density rela- 
tion 



J (Pbulk) — 



Pbulk (1 — Pbulk) 



£- 



l> a Ik 



(11) 



(for 7=1). Exploiting this relation and regarding our 
model as two or three TASEPs joined by slow sites, the 
simplest level of mean- field theories can be built. Ours is 
similar to, but simpler than, the approach in (33| for the 
single defect case. The main difference lies in the match- 
ing condition, i.e., what approximate expression for the 
current across the slow site to use. After comparing the 
two approaches, we proceed to build the case for TASEP 
with two defects. 



A. One slow site 

When a single slow site (q < 1) is located at k, the 
system can be treated as two sublattices: [1,/c] and 
[k + 1,N], referred to as the left and right sublattices, 
respectively. Associated quantities will appear with sub- 
scripts L and R. The two sections are coupled through 
the slow site by having the same current in the steady 
state. Given this constraint, there are only two viable 
scenarios for the sublattices, out of the 3x3 logically pos- 
sible ones: H/L and M/M. 

First, let us consider the H/L case which was one stud- 
ied extensively in (33|. The current for each sublattice 
can be written as: 



Jl 



1 + 13 J 



Jr 



Qfl(l - an) 
1 + a R £ 



(12) 



where (3l and an are the effective exit and entry rates, 
to be determined later. By definition, the entire system 
reaches steady state when = Jr, which yields /?£, = 
an- Of course, these are intimately related to the (bulk) 
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densities through pl — 1 — [3l and pr = ian/ (l + a/jl), 
so that 



i + (i- Pl )£ 



PR 

£ 



Another way to regard this relation is that both densities 
lead to the same current, which we denote by J (a value 
to be determined, and equal to Jl = Jr)- So, the high 
and low densities can be written as p+ (J) and p_ (J), 
respectively, being the two roots to Eqn. (fTTj) . They will 
play a crucial role when we impose the matching condi- 
tion, thereby fixing all quantities as a function of q. 
The exact equation for "matching" is 



J = q{rk(l- n k +Jt)) 



(13) 



in which r k and n k +£ lie in L and R, respectively. Now, 
the right can be expressed as, again exactly, p (k\k + £), 
the probability for finding a ribosome at k, conditioned 
on the presence of a hole at k + £. 

Since we have "broken" the system into two separate 
TASEPs, a naive approximation is to begin with 

Jnmf = q (r k ) ((1 - n k+ t)) 

where the subscript stands for "naive mean field." Re- 
garding this as Eqn. |J7J| for the R sublattice, we have 
a R = QPk- Now, p k is in the L sublattice, and must 
be related to pl in a mean-field approach. The most 
naive assumption is that p\ is the same as its average in 
the bulk, i.e., p r bu i k , which would be Pl/£ in this case. 
However, this turns out to underestimate p (k\k + £) se- 
riously. Indeed, the "pile-up" near a blockage (e.g., in 
Fig. O shows that p\ is significantly higher than its bulk 
value as well as the densities on the £ sites before. Thus, 
we propose that a better approximation would be to re- 
place p\ by pl, and we write 



an = qpL 



(14) 



Using p L = 1 
q/ (1 + q) and 



Pl and f3 L = a R , so that a R = (3 L 



PL = 1/(1 + q), p R = q£/(l + q£), 

we arrive at Jnmf = q/ [(1 + q)(l + q£)]- The premise 
behind this line of arguments is that the system is in 
H/L, so that both aR and (3l should be less than \ ■ 
Therefore, this expression for the current should be valid 
only if it is less than the maximal value (x 2 )- In other 
words, the domain of its validity is limited to q < l/\fi. 
For higher g, this approach predicts that the system will 
be in an M/M phase, with maximal current. Note that 
such a phase cannot occur with a slow defect in the £ = 1 
case, where M/M can be accessed only with q > 1. In 
an earlier study [Hj], the parameters chosen (q = 0.2 
and £ = 12) also precluded the presence of this phase, 
although we believe (see below) that this phase cannot 
be present if the blockage is in the center (A: = N/2) or 



deep in the bulk. We summarize this "naive mean field" 
by 



An alternative approximation for eqn. CE 
posed earlier (33|: 



J = q, 



eff 



PL 



1 



PR \ 



£ J \l- PR £/£) 



was pro- 



(16) 



The last two factors can be recognized as (r k ) and the 
MG approximation for the effective hole density Q ■ The 
first factor is a little more subtle [33|: Considering that 
the transit time for a single particle through the slow 
site (in the absence of steric hindrance) is g -1 +£, q e ff is 
defined as the average rate to move just one step in this 
process: 



Qeff 



q£ 



l + q£ 



The end result for the current is the solution to the alge- 
braic equation 



J = q e ff 



P+ (J) [1 - g_ (J)] 

l-p- (J) £ 



Here, we give an explicit form (which displays the £ = 1 
limit well) 



Jskl — 



Q 



[1-Q + s/1 - 2Q) £ 



(17) 



where 



2g£(l + q£) 



(1 + q + 2q£) 

Note that, for any q < 1, this approach predicts that the 
current is less than the maximal value of x 2 and so, the 
system is always in the H/L phase. 

The results of both mean-field predictions for J as 
a function of q are shown in Fig. 1141 along with two 
sets of data: J q (1) and J q (363). As expected, Jskl 
was purpose-built for two infinite TASEP's connected by 
a slow site and provides a better fit to the data with 
the blockage deep in the bulk (k = 363 in N = 1000). 
On the other hand, it is understandable that, e.g., for 
k = 1, aR must be very close to q. Thus, we may expect 
that the system will have maximal current for q > l/\f£. 
This behavior is confirmed by the data, as illustrated in 
the figure for £ = 12. Since Jnmf (q) has the property 
that it saturates at x 2 for q > x, it provides a better 
fit for J q (1). Of course, we recognize that, as mean- 
field theories, neither (fT4"| nor (|16p are the first step in a 
systematic expansion, so that they may better be thought 
of as "semi-phenomenological" . 
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FIG. 14: (Color online) Comparisons of the current, J, as a 
function of q. The legend labels the two sets of simulation 
data (slow site at k = 1 and 363) and predictions from two 
mean-field approximations. 



More importantly, there is a more serious, inherent 
limitation in this level of mean-field theory. It cannot 
account for the full fc-dependence in J q (k), since it deals 
only with infinite systems. One possibility to incorporate 
finite-size effects (of the left sublattice) in this kind of 
theory is to exploit the MG expression Q 



Jmg 



PiPi+i 



P\+i 



for sites near the boundaries. Using it in a recur- 
sion relation for finding both the density profile and 
the current, we will arrive at a /c-dependent expression, 
Jmg = 1) 0L> k), which replaces Jl in the first expres- 
sion of Eqn. (H2J). However, the high density fixed point 
of this recursion relation is unstable and very careful nu- 
merical analysis will be necessary. Preliminary work on 
this approach is promising and will be reported elsewhere 

We end this subsection by noting that the effects of 
a single defect in TASEP have also been investigated in 
[32| . Unlike our focus here - the dependence of J on 
the location of the slow site, they are concerned with a 
"multicritical system," i.e., a = f3 = 1/2 for the £ = 1 
case. Putting the defect at the center of the lattice, they 
explored density profiles in detail, finding power law tails 
on both sides of the defect with ^-dependent exponents. 
By contrast, our choice of a = = 1 places us far from 
the multicritical point. We have no reason to expect 
similar power laws. 



B. Two slow sites 

The most general TASEP with just two slow sites can 
be quite involved, since the parameter space is four- 
dimensional: {qi, <72, k%, k 2 }. To carry out a manageable 
investigation, we let both sites be deep in the bulk, so 
that only the distance between them, d = (k 2 — k\), plays 
a significant role. Further, as pointed out above, the cen- 
tral section resembles an ordinary TASEP with a and (3 
controlled by q\ and 52, respectively. Therefore, it is the 
smaller (slower) of the two rates which limits that cur- 
rent, which in turn dictates the current through the whole 
system. Thus, we will focus only on the q\ = q 2 = q case. 
Our parameter space will then resemble the single slow 
site case. 

Following the single defect case, the simplest levels of 
mean-field theory treat our system as three subsections 
with obvious labels: L, C, and R. From our discussion, 
only two (out of the many logical possibilities) combina- 
tions of phases, H/S/L and M/M/M, are expected to be 
viable. In addition to Eqn. (fT!?)) . we have 



Jc = 



«c(l - c<c) 
1 + a c l 



(18) 



Since the defect rates are identical, we fully expect that, 
for such a mean-field theory, (3c — etc- Now, matching 
the currents of the subsections, we immediately arrive at 
Jl = Jc = Jr and so, (3 L = a c = Pc = ct R . From 
here, the "naive" mean- field approach for H/S/L pro- 
ceeds identically to the above. The argument relies on 
the presence of a shock in the central section, so that 
there is a low (high) density region near site k\ + 1 (k 2 ) 
and we can impose the same discontinuity in the den- 
sities across both defects, i.e., p + (J) before and p- (J) 
after. Thus, we again arrive at Jnmf (q), given explicitly 
in Eqn. (fT5|) . The same argument can be applied to the 
next level of a mean-field approximation, which predicts 
Jskl (q), as in Eqn. (fT7|) . The major difference between 
the two approaches, as in the single slow site case, is the 
absence of the M/M/M phase in the latter. Meanwhile, 
their limitations are similar: The d dependence in J q (d) 
cannot be accommodated without serious modifications. 

Nevertheless, the spirit of these approximations can be 
exploited to provide J q (1) in the q — > limit. Since the 
central section consists of just one site, there can be no 
shock. Instead, the remnant of the shock is reflected in 
the average density there. It is more convenient to re- 
gard the system as two infinite TASEP 's, with nontrivial 
matching across a "doubly - slow site." Of course, we 
cannot expect to find any of the fascinating profile de- 
tails (e.g., inset of Fig. [TO]) , but we should be able to 
obtain the "coarser" information, such as the currents. 
The goal is to understand the behavior of ^{q — > 0) 
(cf. Fig. I12[) in, say, the first two non-vanishing orders 
in q. 

Now, for q <C 1, we are naturally in the H/L phase and 
the crudest approximation should suffice for the lowest 
order in the current. So, we let the bulk densities be 
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at their extremes (i.e., one and zero) and simply consider 
the time it takes for a particle to move through the block- 
age, from the moment its predecessor is "released." The 
current is just the inverse of this quantity, i.e., 



-+1-2 


-l 

_^ q ' 




~* 2 . 



l-l(£-2) + 



(19) 



where we have included O (g 2 ) terms for computing the 
next order. But, at this order, we should also take into 
account that, occasionally, the density before/after the 
blockage deviates from unity/zero by virtue of the right 
hand side of Eqn. (fTTjl being non-zero. Thus, these den- 
sities are 

p L -> 1- J =l-q/2+ ... 
p R -» J£nq£/2 + ... 

and further suppress the current at the next-to-lowest 
order through the factor 



Combining these factors, we arrive at 



(20) 



If we use exactly the same arguments for the q — > limit 
current in the case of one slow site, we find, instead of 
(USD, 



1+1. 



and, instead of J20 



q[l-q(£-l) + ...} , (21) 



P£(l-p«)->l-g(* + l) + ... • 

Finally since J q [d — * oo) is the same as the single- 
blockage current, we write 



Jq^o (d -> oo) q [1 - 2g£ + 



(22) 



so that 



A 2 (<z^0) 



It is remarkable how well this crude approximation 
agrees with the data in Fig. [T2J There is no doubt 
that all curves extrapolate to the ^-independent value 
of 1/2 at q = 0. As for the slope at the origin, we can 
obtain a good estimate from the lowest q data points, 
using [A 2 (q = 0.02) - 0.5] /0.02. The values obtained 
from simulations for I — 1,2,4,6, and 12 are 0.92, 
1.55,2.53,3.44, and 6.06, respectively. 

We are aware that the expansion (|2"2")) differs from the 
small q limit of Jnmf- Unfortunately it is difficult to 



implement the same scheme for Jnmf here, since we 
must start from the exact pair of equations: 



J = q (r k (1 - n k+e )) = q (r k+ i (1 - n k+ i +1 )) 



(23) 



Various attempts at approximating p k or p k+1 led to 
poorer results. Alternatively, we could exploit the ar- 
gument in SKL [33[ | and consider the average time to 
traverse both slow sites, 2/q + (£—2). This gives us a 
new effective q: 



leff 



q£ 



2 + q(£-2) 



which can be inserted into Eqn. (|16p . The result is 
A 2 (q -> 0) -> § + f (£ + 2) + the O (q) term of which 
differs from the data by about a factor of 2. Clearly, 
mean-field approaches are far from ideal for finding quan- 
titative predictions of J q (d). On the other hand, either 
Jnmf (9) or Jskl (9) provide tolerable results when the 
blockages are from from each other or the boundaries. 
Such variations in the quality of mean-field theories point 
to the importance of correlations. Considerable efforts 
appear to be necessary for a comprehensive, yet relatively 
simple, theory. 

Let us end this section with another method which 
could possibly improve the theoretical predictions [Ell ], 
especially for the first few values of k. The idea is to 
find exact results, by solving the full master equation, for 
TASEPs (with extended objects) on very small lattices 
and then to match these to an infinite system (the R 
sublattice). One expectation is that, as in the £ = 1 case, 
the finite-size current is larger than its counterpart for the 
infinite TASEP at the same (a, /3). This approach may 
eventually provide the essential argument to understand 
the increase in J q (k) as k becomes smaller. Similarly, 
this idea can be applied to the case with two slow sites 
when d is 0(1). Work is in progress to explore these 
avenues. 



SUMMARY AND CONCLUSIONS 



In this study, we consider an inhomogeneous TASEP 
with open boundaries and populated with particles of fi- 
nite extent, £. The hopping rates are uniform (set at 
unity) except for one or two sites ( "defect bonds" ) , where 
the rates, q, are different (faster or slower). We are inter- 
ested in the effects of these local defects on the density 
profiles and the currents through the system. Simulations 
with various £ < 12 show that fast sites have no effect 
on the current, but induce discontinuities in the density 
profiles. In contrast, slow sites generate a significant re- 
duction of the current as well as nontrivial structures in 
the profiles, e.g., long tails behind the blockage, with pe- 
riod £. These findings are entirely consistent with similar 
studies in the past, most of which were restricted to £ = 1 
[3H HH . If the inhomogeneities are deep in the bulk and 
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far from each other, the current depends only on q and 
can be understood through simple mean-field considera- 
tions. Through the current-density relationship (for an 
infinite homogeneous TASEP), the overall densities in 
each of the defect-free sections can be also predicted, so 
that the various "phases" of these subsystems can be un- 
derstood. The distinguishing feature in our study is how 
the location of the defects affects the behavior of the sys- 
tem. For the case of one slow site, the current is slightly 
but measurably enhanced when the defect approaches 
the boundary. On the other hand, a drastic reduction of 
the current is observed when two slow sites are brought 
closer to each other. It is tempting to interpret these ef- 
fects as "interactions" between the defects and to seek a 
formulation that can describe them quantitatively. 

At present, neither the enhancement nor the suppres- 
sion of the currents can be understood in terms of simple 
mean-field theories. The essential limitation is that they 
are based on matching homogeneous TASEPs of infinite 
length. More sophisticated versions, relying on recursion 
relations for the particle density at each site, may be ex- 
ploited to deal with the finite subsections and provide 
some promise for a better understanding of such effects. 
For very small subsections, such as k, d < 5, it may be 
possible to find exact solutions (even for £ > 1) that 
can be used to match the mean-field descriptions for the 
macroscopic subsection(s). Work is in progress to inves- 
tigate these approaches systematically. 

Beyond one or two blockages, we should study systems 
with multiple slow sites, as our eventual goal is to under- 
stand the properties of fully inhomogeneous TASEPs. In 
particular, though restricted to just one or two slow sites, 
our findings - that the relative locations of blockages are 
important - will have implications for translation. For ex- 
ample, they are directly applicable to "designer genes" , 
which consist of many repeats of the same codon, except 
at one or two locations. Using the abundance of associ- 
ated aa-tRNAs as a control for the elongation rate across 



any particular codon, we can test our results directly on 
such genes. Thus, it will be interesting to see the physi- 
cal manifestation of, e.g., enhancement and suppression 
of production rates of such artificial proteins, depend- 
ing on the placement of the slow codon(s). Similarly, 
reproducing the intriguing ribosome density profiles will 
be revealing. More important than "designer genes," we 
should consider the implications for real genes. Our re- 
sults should provide, at the least, some simple qualitative 
insights. We can obviously maximize the production rate 
of a particular protein associated with a certain real gene 
by systematically replacing all slow codons with synony- 
mous, faster ones. However, for most genes this opera- 
tion will require a large number of replacements. Instead, 
with our findings in mind, we can achieve considerable 
increases in the production rates by making only a few 
substitutions, namely, by replacing the slowest codons, 
or a cluster of nearby slow codons. The ratio of cur- 
rent enhancement to the number of codon replacements 
may be used to quantify how "optimal" a certain set of 
substitutions is. This idea can be applied to finding op- 
timal means to suppress prodcution rates as well. Simu- 
lation work with TASEPs associated with real genes is in 
progress and we hope to demonstrate that these concepts 
are viable. 
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I. INTRODUCTION 



A better understanding of non-equilibrium steady 
states in interacting complex systems forms a critical 
goal of much current research in statistical physics. In 
this pursuit, the totally asymmetric simple exclusion pro- 
cess (TASEP) [1-6] has played a paradigmatic role. It 
provides a nontrivial, yet exactly solvable, example of 
phase transitions far from equilibrium, taking place even 
in one-dimensional (ID) lattices. At the same time, it 
also serves as the starting point for the modeling of many 
physical (driven diffusive) processes, such as translation 
in protein synthesis [7-9], inhomogeneous growth pro- 
cesses (e.g. Kardar-Parisi-Zhang growth) [10, 11] and ve- 
hicular traffic [12, 13]. 

In its simplest version, the TASEP involves a single 
species of particles hopping to nearest-neighbor sites, in 
one direction only, along a homogeneous ID lattice. Pro- 
vided the destination site is empty, the rate for the par- 
ticle hop is fixed at 7 (typically chosen as unity without 
loss of generality). With periodic boundary conditions, 
the steady-state distribution is trivial [14] but the full 
dynamics is quite complex [15-17]. With open boundary 
conditions, particles are injected with rate a (in units of 
7) at one end and drained with rate [3 at the other end. 
The competition of injection, transport and drainage in- 
duces a nontrivial phase diagram in the a-(3 plane [1-6], 
reflecting a highly nontrivial steady state. Three phases 
are present: a maximum-current phase for ev.,/3 > 1/2, 
and a low- (high-) density phase for a < (3, a < 1/2 
{(3 < a, (3 < 1/2). Not surprisingly, there are also rich 
dynamical aspects [18, 19]. 

To model protein synthesis, each site on the lattice 
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represents a codon on the messenger RNA (mRNA), and 
the particles represent the ribosomes. Injection, hop- 
ping, and drainage are associated respectively with ini- 
tiation, elongation, and termination in biological terms. 
The quantity of interest, namely, the (steady-state) pro- 
tein production rate, is identical to the (stationary) par- 
ticle current. Clearly, the simple TASEP falls short of 
the biological system in several significant aspects. One 
is that an individual ribosome "covers" several codons 
[7, 20, 21], as opposed to a particle occupying only a 
single site. Another is that, in all naturally occurring 
mRNAs, the codons carry genetic information and there- 
fore necessarily form an inhomogeneous sequence. Thus, 
the elongation rate of a ribosome is unlikely to be uni- 
form; instead, the hopping rate, ji, of a particle becomes 
a function of the site i. For example, it is well known 
that translation slows down at specific codons (see, e.g. 
[9, 22-25]), with potentially significant consequences for 
protein production rates. Indeed, the steady-state cur- 
rent may depend sensitively on not only the frequency of 
each codon's occurrence, but also the order of their ap- 
pearance in the sequence. Both of these issues - extended 
objects and inhomogeneous rates - have been addressed 
recently in separate contexts which we summarize briefly 
in the following. 

The results associated with inhomogeneous (quenched 
random) rates fall into two broad categories, in the sense 
that the randomness can be associated with the parti- 
cles [26-28] or with the sites. Randomness of the former 
type is more relevant for vehicular traffic where it ac- 
counts for a variety of driver preferences. In contrast, 
the disorder in the protein case is clearly site-dependent, 
leading to spatially non- uniform hopping rates 7$. Re- 
stricting ourselves to this class, we can consider the ef- 
fect of having a whole distribution, or very specific con- 
figurations, of 7^. Starting from given distributions, two 
groups [29, 30] studied the resulting disorder-average in 
periodic systems. To mention just one significant effect, 
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the current-density diagram develops a plateau: limited 
by the smallest rate in the system, the current becomes 
independent of density over a range of densities. Harris 
and Stinchcombe [30] also extended this work to open 
systems. While these studies may be of some interest to 
mixtures of many different mRNAs, our primary inter- 
est here is to understand how the production rate of a 
specific protein is associated with a specific genetic se- 
quence. As a first step towards a solution, we adopt the 
approach of several other studies [9, 31-33], by focusing 
on the effects of a few localized inhomogeneities, i.e., hop- 
ping rates which are uniform except at a handful of sites 
[34]. As a synthesis of these studies, we will explore in 
some detail the consequences of having extended objects 
and locating one or two slow sites at a variety of positions 
on the lattice. In this manner, by introducing more and 
more sites with a range of rates, we hope to understand 
inhomogeneities in a systematic way, setting the stage for 
further investigations of the translation process. 

A full comprehension of the effects of slow sites on 
the particle current may have potentially significant ap- 
plications in biotechnology. While there are 64 distinct 
codons, proteins are chains composed of just 20 amino 
acids. So, many different mRNAs (codon sequences) 
can code for the same particular protein. Moreover, the 
amino acid is incorporated into the growing chain by 
an important intermediary, the so-called transfer RNA 
(tRNA), which carries the complementary anticodon. It 
turns out that the mapping between codons and tR- 
NAs is also not precisely 1-1. For example, in E. coli, 
the genetic code actually involves 61 sense codons and 
about 46 tRNAs with associated anticodons [35]. Mean- 
while, for a given mRNA sequence, the protein produc- 
tion rate is often modeled in terms of (generally accepted) 
charged-tRNA (aminoacyl-tRNA, or aa-tRNA) concen- 
trations [24], so that different sequences can result in 
different production rates for the same protein. By elu- 
cidating how the spatial distribution of defects, especially 
of bottlenecks, affects translation rates, we can pinpoint 
those clusters of codons which are likely to have the most 
significant effect on the production rate of the associated 
protein. Exploiting the degeneracy in the mapping from 
mRNA sequence to protein, we can provide guidance as 
to how a few selected, local modifications of the mRNA 
can optimize the production rate of a given protein. 

Our paper is organized as follows. In Section II, we 
define the model and provide a more detailed description 
of previous work, concerning exclusion processes with ex- 
tended objects or spatially inhomogeneous rates. A brief 
discussion of a previous mean-field analysis is also in- 
cluded. In Section III, we present our Monte Carlo re- 
sults. We focus especially on the implications of hav- 
ing extended objects by varying the particle size. We 
first consider the interaction between one slow site and 
the system boundary, and, motivated by the resulting 
findings, turn to the interactions between two slow sites. 
This provides new insights for genes containing clusters 
of slow codons, which occur frequently in, e.g., E. coli, 



Drosophila, yeast and primates [9, 36, 37]. In Section IV, 
a complete investigation of systems with inhomogeneities 
is presented using a mean-field approach. Section V con- 
tains our conclusions and a summary of open questions. 

II. MODEL SPECIFICATIONS AND KNOWN 
RESULTS 

The TASEP is defined on a ID lattice of N sites. We 
introduce an index i = 1, 2, N to label the sites. Each 
site (codon) is either occupied by a single particle (ribo- 
some) of length I (in units of sites) or empty. A micro- 
scopic configuration of the system can be uniquely char- 
acterized in terms of a set of occupation variables, \n{\, 
taking the value 1(0) if site i is occupied (empty). Of 
course, the extended nature of the particle induces strong 
correlations in {rii}, in the sense that a single ribosome 
always covers I consecutive sites. Yet, at any given time, 
only one of the covered codons is being "read" (i.e., the 
codon is "covered by the aminoacyl site", or A site, of 
the ribosome) and translated into an amino acid. Here, 
we refer to the associated location on the ribosome as 
the "reader" (of the genetic code). For our purposes, it 
is not essential which one of the £ sites is labeled as the 
reader, and so we follow the convention in [8] and choose 
the first (leftmost) site. Hence, the statement "a ribo- 
some (or particle) is located at site i" implies that the 
reader is located at site i and the subsequent 

l = l-\ (1) 

sites are also occupied. Naturally, the position of the 
reader determines the elongation rate, i.e., ji, since the 
ribosome must wait for the arrival of the aa-tRNA with 
the i-specific anticodon before it can move to the next 
site. Clearly, the reader locations can also be used to 
label a microscopic configuration, i.e., we can define the 
reader occupation number at site i as r,. The sets {n^} 
and {ri} are uniquely related to each other. Moreover, 
due to the extended size of a particle, strict constraints 
are built in (e.g., = 1 implies r i+ i = ... = r i+ £-i = 
and Hi = ... = n i+ e-i = 1). As a consequence, neither 
set can be arbitrary and serious correlations arise as soon 
as i > 1 [38]. 

In our simulations, we adopt a random sequential up- 
dating scheme and keep a list of locations of readers. In 
addition, the site i = is always occupied by a "virtual 
reader," which accounts for particles entering the system 
(initiation). At the beginning of each Monte Carlo step 
(MCS), we first find the number of particles in the sys- 
tem and label it M. Then, we randomly select an entry 
from this list of M + 1 readers. If the chosen reader is 
virtual (i.e., i = 0), a new particle enters the lattice with 
probability a, provided all the first i sites are empty. If 
the chosen reader is real, say, at site i > 0, the associ- 
ated particle is then moved to site i + 1 with probability 
7i, provided site i + t is empty. With this notation, we 
can also write the initiation and termination probabilities 



3 



(a) 



a 



z^msssm 

|- k H 



P 



(b) 



a A A 
, 1 1 , i i^s^m 



v- d h 



| | | | 



FIG. 1: Sketch of a TASEP for particle size ^ = 6 with (a) a 
single slow site at position k, with rate q, and (b) two slow 
sites with rate q, separated by a distance d. 



(a and (3) as 70 and 7 at, respectively. To be complete, 
the sites beyond the lattice are by definition empty, so 
that once a particle reaches N — £ + 1, it will not ex- 
perience steric hindrance (see Fig. 1 for a sketch of this 
process). These processes have been termed "complete 
entry" and "incremental exit" [39] . Other entry and exit 
rules can be considered, but are believed to be incon- 
sequential provided £/N <C 1. Each MCS consists of 
M + 1 such attempts, giving an even chance, on average, 
for each particle (ribosome) in the system to elongate or 
terminate, as well as for an initiation event to occur. 

Starting with an empty lattice, we typically discard 
2 x 10 6 MCS to ensure that the system has reached the 
steady state. Unless otherwise noted, good statistics re- 
sult if we average over least 2 x 10 4 measurements, sep- 
arated by 100 MCS in order to avoid temporal corre- 
lations. Such steady state averages will be denoted by 
(...). To reduce the number of parameters in the model, 
we study systems with a = (3 = 7$ = 1, except at one or 
two sites. The system sizes (N) range from 200 to 1000, 
with most data taken from TV = 1000. 

To characterize the state of the system, we monitor 
several observablcs. The most obvious is p\ = (rj), a 
quantity we will refer to as the ribosome (or reader, or 
particle) density. Of course, J2i Pi 1S J ust tne average 
number of particles in the system (i.e., ribosomes on the 
mRNA). Thus, the overall particle density N'^-J^iPi 
is bounded above by l/£. Another interesting variable 
Pi = (n>i), labeled as the "coverage density", is the prob- 
ability that site i is covered by a particle (regardless of 
the location of the reader). Needless to say, the pro- 
file for the vacancies is given by the local hole density, 
p\ = 1 — pi. The overall coverage density, A^ 1 p,, 
may reach unity and provides a good indication of how 
packed the system is. The two profiles are related by 
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A quantity of great importance to a biological system 
is the steady-state level of a given protein. If wc assume 
that the degradation rates are (approximately) constant, 
i.e., independent of protein concentration, then these lev- 
els are directly related to the protein production rates. In 
our model, such a rate is just the average particle current 
J, defined as the average number of particles exiting the 
system per unit time. In the steady state, it is also the 
current measured across any section of the lattice. For 
simplicity and to ensure the best statistics, we count the 
total number of particles which enter the lattice over the 
entire measurement period (at least 2 x 10 6 MCS in most 
cases). 

In this study, we focus on two simple types of inhomo- 
geneities: one or two "slow" sites (Fig. 1). Their locations 
specify the only inhomogencitics in the rates. 

One slow site, at position k. We denote jk by q (< 1). 
This corresponds to a bottleneck in the lattice. We are 
especially interested in the dependence of the current, 
denoted by J q (k), on the parameters q and k. 

Two slow sites, at positions k\ and ki with separation 
d= (k 2 — fci) and rates q\ t 2 = 7fei, 2 - We find that, when 
qi qi > the current is controlled mainly by the smaller of 
the two, with little dependence on d, in agreement with 
a simple mean-field theory to be discussed in Section IV. 
Therefore, most of our attention will be devoted to the 
case with q\ = qi = q < 1. Moreover, we choose to limit 
our study to both sites being far from the boundaries. 
Then, the current is insensitive to their average position 
(fc 2 +fci)/2, and we can investigate J q (d). Note that these 
are precisely the systems studied in [9], except that we 
consider particles with a range of sizes: I — 1,2, 4, 6, and 
12. While there are qualitative similarities, we will dis- 
cuss the quantitative differences due to I > 1, as well as 
interesting phenomena associated with the density pro- 
files. 

Let us now provide the context of our work by briefly 
reviewing some related earlier studies. The homoge- 
neous case (71 = ... = 7at_i = 1) with £ = 1 is ex- 
actly soluble [2-6], and displays three phases in the a-(3 
phase diagram. For £ > 1, no exact solutions exist [40]. 
Analytic approximations using various mean-field ap- 
proaches [7, 8, 20, 39, 41] predict the presence of the 
same phases, though the phase boundaries depend on £ 
(Fig. 2) through the combination [8] 
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with the understanding p\ = for i < 0. 



Monte Carlo studies [8, 39] largely confirm these conclu- 
sions. 

The three phases carry different currents and dis- 
play distinct density profiles [2-6, 8, 33, 39]. Apart 
from "tails" near the boundaries, the (coverage) den- 
sity profiles approach uniform bulk values in the ther- 
modynamic limit, i.e., pi — > Pbuik, for 1 <C i <C N. 
For a < x an d a < [3, the system is in a low-density 
phase (L), characterized by pbuik — £ct/{l + al) and 
J = a(l — a)/ (l + at). A high-density phase (H) pre- 
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FIG. 2: Phase diagram for an ordinary TASEP. On the dashed 
line, the H and L phases coexist. 



TABLE I: J-pbuik relation for particles of size £ (£ : 
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phase 



current J 



bulk density p bu i k 



L 
H 

M 



a(l-a)/ (1 + at) 
13(1-/3)/ 



£a /(l + at 
1-/3 

i-x 



vails for [3 < x and (3 < a, with bulk density pbuik = 1 — /? 
and current J = /3(1 - /3)/ (l + /3l). For > x, 

the system is in a maximum-current phase (M), where 
Pbuik = 1 - X and J = x 2 - On the a = (3 < \ line 
(dashed line in Fig. 2), the system consists of two macro- 
scopic regions, characterized by a low (high) density re- 
gion near the entry (exit) point. The two regions are 
joined by a shock front that performs a random walk. 
This is often referred to as the "shock phase" (S). Ta- 
ble I summarizes the J — pbuik relation for TASEP with 
extended objects. There is good agreement between sim- 
ulations (with I < 12) and analytic results for these bulk 
quantities [8, 39]. The details of the profile for £ > 1, 
especially near the lattice boundaries, are less well un- 
derstood. While periodic structures (of period £) can 
be expected, mean-field theories [7, 20, 39] were success- 
ful in capturing only a limited part of the phenomena 
observed. We will return to these considerations in Sec- 
tion IV. Beyond homogeneous systems, several studies 
introduced one or more "impurities" into TASEP with 
periodic boundary condition. A single "slow" site in- 
duces a shock in the density profile with some interesting 
statistics [42-45, 47]. Subsequently, generalizations to 
systems with a finite fraction of slow sites, randomly lo- 



cated, were also investigated [29]. For the richer case of 
the open boundary TASEP [9, 31-33, 48], Kolomeisky 
focused on point particles (£ = 1), with a single impu- 
rity at the center of the lattice [31], so as to mimic a 
defect situated deep in an infinitely long system. The 
consequences of the defect having both faster (q > 1) 
and slower (q < 1) rates were explored. By matching 
two ordinary TASEPs across the defect, the properties 
of such systems in the a-(3 plane can be well described 
[31, 48]. While a fast site has no effect on the phase 
diagram, a slow site leads to a shift of the M-H and M- 
L phase boundaries to ^-dependent, smaller values of a 
and f3. The density profiles are quite sensitive to the ex- 
istence of a defect site [9, 31]. Kolomcisky's approach 
was generalized to the I = 12 case in [33], with similar 
levels of success. Below, we will provide further details 
of this work, on which we base much of the analysis of 
our problem. Ha and den Nijs also studied the i = 1 
open boundary TASEP with a single defect at the center 
[32]. Focusing on the multicritical point a = (3 = 1/2, 
they were mainly interested in the so-called "queueing 
transition" and its critical properties. Detailed results of 
density profiles, such as power law behavior and critical 
exponents, were obtained in the region q = q c . Here, 
q c denotes the critical value of q below which the bulk 
density in front of the slow site deviates from the den- 
sity behind the blockage. By contrast, our focus here is 
essentially that of [9, 48], namely, how does the number 
and the locations or spacings of the slow sites affect the 
current through the system? In the single slow site case, 
[9] investigated mainly the overall current as a function 
of q, while [48] analyzed the current as a function of k 
(the distance between the defect and the entry point). 
For the case of two defects, both studies find that the 
spacing between them plays a significant role for the cur- 
rent. A finite-segment mean- field theory in [9] provides 
excellent agreement with data. In particular, clustered 
defects reduce the current much more effectively than 
well separated ones. In this sense, we can regard these in- 
vestigations as exploring the "interactions" between the 
slow site(s) and/or the boundaries. Since both studies 
are restricted to point particles (£ = 1), our intent is 
to explore the effects of having extended objects (with 
£ < 12). Though we expect qualitatively similar behav- 
ior, as pointed out in [9], we also find noteworthy quanti- 
tative differences. Finally, we should mention that larger 
numbers of defects do not lead to significantly different 
effects [9], so that we limit ourselves to one or two slow 
sites here. 



III. MONTE CARLO RESULTS 

In this section, we present our Monte Carlo results. 
For convenience, we use a consistent color coding scheme 
(online only) for the various particle sizes, as specified in 
Table II. The data here consist of the overall currents 
and the density profiles of both coverage and ribosomes. 
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TABLE II: Color coding scheme 



size £ 


online color 


1 


black 


2 


red 


4 


brown 


6 


green 


12 


blue 



Our focus will be how these quantities depend on q, k (for 
the case with a single defect), and d (for the case with 
two defects) . Although the profiles arc difficult to extract 
experimentally, the reader profiles will be of interest in 
subsequent studies involving real gene sequences, since 
they provide information on how frequently a particular 
tRNA is bound to the mRNA. By contrast, the currents 
are easily measurable and our results here may generate 
more immediate interest. 



A. One slow site 

We begin by placing one slow site on the lattice as in 
Fig. 1(a). Fig. 3 shows several coverage density profiles 
for a typical choice of parameters: N = 1000, q = 0.2, 
and k = 82 with £ = 1,6, 12. As expected, we observe 
pile-ups of particles due to the blockage - a high (low) 
density region before (after) the bottleneck in all three 
cases. However, due to the lack of ordinary particle-hole 
symmetry in the £ > 1 cases, the average densities on 
either side of the slow site are not symmetric around 
0.5. Instead, they are roughly related through the J-pbuik 
relations in the H and L phases, summarized in Table I. 
In detail, the profiles are quite different: The "tails" , i.e., 
the deviations from the bulk values, are quite noticeable 
in the vicinity of both the slow site and the edges of the 
system for the £ — 6,12 cases. The inset exposes more 
clearly that there are period £ structures in the profiles 
[7, 20, 39], especially just before the slow site. 

A more dramatic difference between point particles and 
extended objects emerges when we plot the ribosome den- 
sity p\, in Fig. 4, corresponding to the inset in Fig. 3. 
Similar to profiles in [7, 20, 39], we find distinct period £ 
structures before the slow site. While the reader "waits" 
to pass the blockage, the readers of the following parti- 
cles tend to catch up and pause at sites k — n£, where 
n = 1, 2, ... The "tails" are even more marked than those 
in Fig. 3. To emphasize the difference between the reader 
and coverage profiles (p\ and p^), we show a case with 
q = 0.05, k = 948, £ = 12 in Fig. 5. Though both profiles 
contain the same information, we see that p\ (lower plot) 
is far more sensitive than pi (upper plot) in showing the 
very long tails (~ 1000 in this example) hidden in the col- 
lective behavior of the particles. At present, the crucial 
ingredients that control the characteristic decay length of 
the pJJ-envelopes have not yet been identified. Certainly, 
these very large length scales are completely absent from 
the £ = 1 systems deep within the H/L phases. 
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FIG. 3: (Color online) Coverage density profiles with one slow 
site of q = 0.2 at k = 82. I = 1,6, 12 (from bottom to top 
in both subsections) and N = 1000. The inset is a magni- 
fied view of the i £ [1, 150] interval, to expose the period £ 
structures. 




FIG. 4: (Color online) Ribosome density profiles with one 
slow site of q = 0.2 at k = 82. t = 1,6,12 (from top to 
bottom in both subsections) and N = 1000. Only the first 
150 lattice sites are shown. 



As for the current, Fig. 6 illustrates its dependence on 
q, fc, and £. Not surprisingly, the current is limited by 
the bottleneck and therefore varies monotonically with q. 
It is also reduced if the particle size increases, an effect 
that can be traced mainly to the particle density being 
effectively lower by the factor £. For point particles, the 
current is not very sensitive to the location of the slow 
site. The enhancement as k approaches the boundary 
of the system - referred to as the "edge effect" [48] - is 
quite small. For larger £, the enhancement is much more 
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FIG. 5: (Color online) Coverage density profile (top) and ri- 
bosome density profile (bottom) with one slow site of q — 0.05 
at k = 948. I = 12 and N = 1000. 
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FIG. 6: (Color online) J q (k) as a function of the location k 
of the slow site for q — 0.2 (lower set of squares); 0.3 (middle 
circles) and 0.4 (upper triangles), (a) I = 1; (b) I — 2; (c) 
I = 6; (d) I = 12. In all cases, JV = 1000. 



pronounced, especially for smaller q. Whatever the mag- 
nitude, in all cases the current increases monotonically 
as the slow site is located closer and closer to the entry 
point. For 1=1, particle- hole symmetry is manifest in 
the microscopic dynamics, so that the symmetry of J q (k) 
under k — ■> N + 1 — k inversion, is obvious [48] . For £ > 1, 
the density profiles confirm the lack of this particle-hole 
symmetry very clearly. Correspondingly, there is a sys- 
tematic asymmetry in the current: J q {k) = J q (N + 1 — k) 
is satisfied only for k < I. The origin of this behavior is 
not well understood. 

The edge effect, and specifically its dependence on q 



FIG. 7: (Color online) Ai(q) for I = 1,2,4,6 and 12. 



and £, can be quantified by the ratio: 



Ai(?) = 



J q (k 



1) 



J q (k — > OO) 



(4) 



Fig. 7 shows that A-i(q) depends on q in a nontrivial 
way. The maxima of Ai(g) occur at lower values of q 
as £ increases, reminiscent of the behavior of the phase 
boundary between M and L/H. With appropriate scaling, 
the curves of A\(q) can be collapsed for large i's. From 
the biological perspective, the edge effect is not easily ob- 
servable since the current enhancement is less than 10% 
for the relevant £. 

Returning to Fig. 6, we note that significant deviations 
from the asymptotic value, J q (po), are found only for fc's 
within a small distance from the boundaries. At a ca- 
sual glance, this range appears to depend on both q and 
£. On closer examination of, say, the most prominent 
case here: (q,t) = (0.2,12), we find that the decay of 
J q (k) into J 9 (oo) fits an exponential quite well (Fig. 8), 
i.e., J q (k) — J g (oo) oc exp (— k/5), with 6 s» 10. Assuming 
this behavior persists in the other cases, we can study the 
(q, £) dependence of this characteristic length and denote 
it by S(q,£). We believe that the origin of this length 
scale can be traced to the presence of tails in the density 
profiles near the lattice boundaries: As the slow site ap- 
proaches the entry point, these tails may seriously affect 
the injection process, and thus, the current as well. For 
the £ = 1 case, we observe that S(q,£) ~ C(^e//)j where 
£ is the characteristic length associated with the bound- 
ary layer of the density profile in the ordinary TASEP. 
Specifically, with entry rate a = 1 and exit rate /?, the 
system is in the H-phase and the profile decays expo- 
nentially into the bulk, as p x — pbuik oc exp(— x/£) with 
[6]: 
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FIG. 8: (Color online) Dependence of J q (k) on k obtained 
from simulation is plotted in squares and the line is a linear 
fit with slope equals -0.11. q = 0.2, £ = 12 and N = 1000. 



Here, the left half of our system is such a TASEP, except 
that we have an effective ft: /3 e // = e?/(l + q) [33]. Us- 
ing these arguments on the three g's shown, we estimate 
decay lengths of about 5 (q = 0.4), 3 (q = 0.3), and 2 
(q = 0.2) lattice constants. Though the data on the dif- 
ferences Jq(k) — Jq(oo) are small and noisy, simulation 
results are consistent with S(q, 1) ~ £.(Peff)- However, 
for £ > 1, there is no analytic result for the boundary 
layers of the density profiles. Moreover, the data sug- 
gest that they are quite complex (e.g., in Figs. 4 and 5). 
Thus, it is unclear how to quantify the picture for point 
particles to the general case of S(q, £). At present, a com- 
plete understanding of both "boundary layers" - in the 
density profiles and in J q (k) - remains elusive. 

If we consider the edge effect as an "interaction" be- 
tween the slow site and the lattice boundaries, the nat- 
ural next step is to explore the interactions between two 
slow sites. In order to avoid edge effects, we place the 
two slow sites sufficiently far away from the boundaries 
and vary their separation. 



B. Two slow sites 

As mentioned in the previous section, the currents in 
the qi q2 cases are essentially controlled by the slower 
of the two rates and so, may be regarded as systems with 
a single slow site. These profiles can be interesting, but 
we choose to restrict our attention here to a study of the 
qx = qi = q case, in which the currents show a nontrivial 
dependence on d, the distance between the two slow sites. 
With two bottlenecks, the system consists of three sec- 
tions: before the first blockage, in between the two, and 
after the second defect. Of course, for small q, the overall 
density in the first (last) section is expected to be high 
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FIG. 9: (Color online) Coverage density profiles for two slow 
sites with q = 0.2. I = 12, d = 100; £ = 6, d = 125; 1 = 2, 
d — 150; and I — 1, d = 170 (Curves are plotted from left to 
right in the mid-subsection and top to bottom elsewhere). In 
all cases, N = 1000. Inset: q = 0.6, I = 12, d = 2400 and 
TV = 7200. 



(low). In these cases, the effective entry and exit rates for 
the central section are also low, so that a wandering shock 
should be present. Hence, the average profile should be 
linear for £ = 1 (and essentially so for larger £ [8]) with 
a positive slope. This behavior is understandable, since 
the section between the two defects is comparable to an 
ordinary TASEP with small a — [3. These expectations 
are generally confirmed by simulations with q < 0.5 and 
various £'s up to 12. Fig. 9 shows typical coverage pro- 
files, for a relatively small rate of q = 0.2. The system 
appears to make a transition from this H/S/L phase to 
an M/M/M phase as q increases. The center profiles be- 
come essentially flat, as illustrated in the inset (where 
q = 0.6 and £ — 12). Details of this transition are being 
explored. 

More interesting are the finer features of the profiles 
in the small q cases. As in the single defect system, the 
profiles exhibit period £ structures near the slow sites. 
To resolve these more clearly, we plot the reader density 
profiles in Fig. 10. In all cases that involve extended par- 
ticles (£ > 1), the readers clearly pile up behind the slow 
sites. Apart from these "jams," another feature emerges, 
namely, a sequence of depletion zones, each of which pre- 
cedes one of the period £ peaks. For £ = 2, the differences 
between the upper and the lower envelope are especially 
dramatic. More remarkably, when the blockages are sep- 
arated by small d's, two different, "overlapping tails" are 
created, as illustrated in the inset of Fig. 10, where d = 1, 
q = 0.2, and £ = 12. Indeed, there are further interesting 
structures for d < £, which will be presented elsewhere. 

Compared to these remarkable characteristics in the 
profiles, the behavior of the currents seems lackluster. In 
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FIG. 10: (Color online) Ribosome density profiles with two 
slow sites of q = 0.2. I = 12, d = 100; t = 6, d = 125; t = 2, 
d = 150; and I = 1, d = 170 (Curves are plotted from left 
to right in the mid-subsection and bottom to top elsewhere). 
Inset, 1 = 2 and d = 1. In all cases, N = 1000. 



Fig. 11, we plot four sets of currents 



J q (d), associ- 
ated with £ = 1,2,6, and 12. In all cases, we see that 
J is considerably suppressed when d is reduced. When 
the slow sites are very far apart, the current behaves as if 
there is only one slow site, consistent with expectations 
from mean-field theories. At the other extreme, when 
the two defect sites are nearest neighbors, the current 
reaches its minimum. Not surprisingly, period £ struc- 
tures emerge as d is varied, illustrated in the inset of 
Fig. 11(d), but become less prominent for d > 50. These 
plots also reveal that, unlike the dependence on k above, 
there are serious deviations from the d — > oo values when 
d is decreased. To quantify this deviation, we define 



A a (g) 



J q (d= 1) 

Jq(d — * OO) 



(6) 



and plot this quantity vs. q in Fig. 12. In contrast to 
Ai(g), we observe that A2(q) exhibits a sizable depen- 
dence on q, especially for small values of q. In the limit 
of q — > the current decreases by a factor of 2! In the 
following section, we will see that this factor can be un- 
derstood via a mean-field approach. 

To summarize our simulation results, two bottlenecks 
near each other have a dramatic effect on the current. We 
may regard this phenomenon as an "interaction" between 
the two slow sites, inducing far more "resistance" when 
they are close than when they are well separated. 

Two additional comments are in order. First, we re- 
turn to one of the predictions of the mean-field theory, 
namely that a second slow site, spaced far apart from its 
partner, should have no further effect on the current. Our 
data indicate that the current for two slow sites, spaced 
far apart, is systematically lower than the current for a 
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FIG. 11: (Color online) J q (d) as a function of the separation d 
between the two slow sites for q = 0.2 (lower set of squares); 
0.3 (middle circles) and 0.4 (upper triangles), (a) I = 1; 
(b) I = 2; (c) I = 6; (d) £ = 12. The inset in (d) is a 
magnified view of the d £ [1, 60] interval, to expose the period 
£ structures. In all cases, N = 1000. 
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FIG. 12: (Color online) A 2 (<?) for £ = 1,2,4,6 and 12. 



single slow site, but only by a very small amount (less 
than 1%). Second, we can again attempt to identify a 
length scale which controls how J q {d) approaches J q (oo), 
as d increases. Since the central section of the system 
displays a shock, it is natural to ask whether the intrin- 
sic width of the shock sets this length scale. According to 
[42, 45], this width covers only a few lattice spacings in 
the periodic TASEP with a single defect. Here, however, 
it appears that the shock is much broader. For example, 
the averaged profile of the shock for the case of q = 0.2 
with point particles is shown in Fig. 13, as well as a sim- 
ple fit using a tanh function [46] with width of about 10. 
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FIG. 13: (Color online) The open circles mark the average 
profile of a shock between sites 149 and 349 with q = 0.2, 
compiled from of a very long run (3 x 10 s MCS) with an 
N — 1000, 1 = 1 system. Details of how raw profiles are 
shifted (so that the shock is located at site x = shown 
here) will be published elsewhere [49]. A simple fit using 
A + B tanh (x/10) is also shown: solid line. 



Intriguingly, this length appears to be comparable to the 
one appearing in Fig. 11(a). More work is needed to fully 
explore these issues. 



IV. MEAN-FIELD THEORETIC APPROACHES 

Mean-field theory is known to agree well with the ex- 
act results for a number of macroscopic quantities in 
the steady state of the I = 1 TASEP (see, for example 
[6]). For extended particles, no exact solution is avail- 
able [40] so that mean-field (and more sophisticated clus- 
ter) approximations form the only route toward some un- 
derstanding of the system's behavior. However, there are 
many levels of "mean- field" approximations [7, 8, 33, 39], 
corresponding to neglecting different types of correla- 
tions. For certain quantities (e.g., currents in large sys- 
tems), predictions from the simplest level are very close 
to the simulation results. For others (e.g., some reader 
profiles), only the most sophisticated level performs ad- 
equately. In all cases, no level of mean-field theory can 
give a good fit to both the current and the profile. A 
thorough discussion of all of these schemes is quite in- 
volved and will be provided elsewhere [49]. Here, we will 
restrict ourselves to the simplest method and compare its 
performance to the simulations. 

All approaches start with the exact expressions for the 
current 

J = a(l-m) (7) 
= H (n (1 - n i+i )) ; i€[l,N-£] (8) 



it [n) 



i € 



[N-e,N-i] 



(9) 
(10) 



In the absence of the steady-state distribution, the most 
naive approximation is to replace (riUj) by (r^) (rij). Un- 
fortunately, the constraints due to particles with I > 1 
are so severe that this approximation is entirely inad- 
equate when j = i + i. Even for the simple case of 
TASEP on a periodic ring, it leads to an erroneous ex- 
pression for J (except if I = 1). Instead, the average 
(coverage) density at site % + I is much larger than the 
(conditional) probability that it is actually covered given 
that the reader is at site i. MacDonald and Gibbs (MG) 
proposed [7] a much better approximation: 



(n (1 - n i+ i)) 



Pii 1 - Pi+e) 



P\p\+i 



1 



p l+ i + p\ +l p\ +l + p\ +t 



As discussed in Section II, the densities far from bound- 
aries are uniform (e.g., pf—Kx, — > Pbulk /(■) and this fact 
provides a good description of the current-density rela- 
tion 



J (Pbulk) — 



Pbulk (1 — Pbulk) 



£- 



l> a Ik 



(11) 



(for 7=1). Exploiting this relation and regarding our 
model as two or three TASEPs joined by slow sites, the 
simplest level of mean- field theories can be built. Ours is 
similar to, but simpler than, the approach in [33] for the 
single defect case. The main difference lies in the match- 
ing condition, i.e., what approximate expression for the 
current across the slow site to use. After comparing the 
two approaches, we proceed to build the case for TASEP 
with two defects. 



A. One slow site 

When a single slow site (q < 1) is located at k, the 
system can be treated as two sublattices: [1,/c] and 
[k + 1,N], referred to as the left and right sublattices, 
respectively. Associated quantities will appear with sub- 
scripts L and R. The two sections are coupled through 
the slow site by having the same current in the steady 
state. Given this constraint, there are only two viable 
scenarios for the sublattices, out of the 3x3 logically pos- 
sible ones: H/L and M/M. 

First, let us consider the H/L case which was one stud- 
ied extensively in [33]. The current for each sublattice 
can be written as: 



Jl 



1 + 13 J 



Jr 



«fl(l - aw) 
l + awl 



(12) 



where (3l and an are the effective exit and entry rates, 
to be determined later. By definition, the entire system 
reaches steady state when — Jr, which yields /?£, = 
an- Of course, these are intimately related to the (bulk) 
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densities through p L = 1 — f3 L and pn — £chr/ (l + (xr£) , 
so that 



1 + 



PR 

£ 



Another way to regard this relation is that both densities 
lead to the same current, which we denote by J (a value 
to be determined, and equal to Jl = Jr). So, the high 
and low densities can be written as p + (J) and p_ (J), 
respectively, being the two roots to Eqn. (11). They will 
play a crucial role when we impose the matching condi- 
tion, thereby fixing all quantities as a function of q. 
The exact equation for "matching" is 



J = q(r k (1 - n k+l )) 



(13) 



in which r k and n k +e lie in L and R, respectively. Now, 
the right can be expressed as, again exactly, p (k\k + £), 
the probability for finding a ribosome at k, conditioned 
on the presence of a hole at k + £. 

Since we have "broken" the system into two separate 
TASEPs, a naive approximation is to begin with 

Jnmf = q {r k ) ((1 - n k+l )} 

where the subscript stands for "naive mean field." Re- 
garding this as Eqn. (7) for the R sublattice, we have 
a R = QP k - Now, p\ is in the L sublattice, and must 
be related to p L in a mean-field approach. The most 
naive assumption is that p\ is the same as its average in 
the bulk, i.e., p\ u i k , which would be pl/£ in this case. 
However, this turns out to underestimate p (k\k + I) se- 
riously. Indeed, the "pile-up" near a blockage (e.g., in 
Fig. 5) shows that p k is significantly higher than its bulk 
value as well as the densities on the I sites before. Thus, 
we propose that a better approximation would be to re- 
place p\ by p L , and we write 



ctR = qpL 



(14) 



Using p L = 1 - (3 L and (3 L = a R , so that a R = (i L = 
q/(l + q) and 

p L = l/(l + g), p R = q£/(l + qe), 

we arrive at Jnmf = q/ [(1 + ?)(!+ q£)\- The premise 
behind this line of arguments is that the system is in 
H/L, so that both a R and /?£ should be less than \ ■ 
Therefore, this expression for the current should be valid 
only if it is less than the maximal value (x 2 ). In other 
words, the domain of its validity is limited to q < 1/ VI. 
For higher q, this approach predicts that the system will 
be in an M/M phase, with maximal current. Note that 
such a phase cannot occur with a slow defect in the £ — 1 
case, where M/M can be accessed only with q > 1. In 
an earlier study [33], the parameters chosen (q = 0.2 
and £ = 12) also precluded the presence of this phase, 
although we believe (see below) that this phase cannot 
be present if the blockage is in the center (k — N/2) or 



deep in the bulk. We summarize this "naive mean field" 

by 

Jnmf - j ±2 for ^ ^ . (15) 

An alternative approximation for eqn. (13) was pro- 
posed earlier [33]: 



<leff 



PR \ 



i - pr£/£J 



(16) 



The last two factors can be recognized as (r k ) and the 
MG approximation for the effective hole density [7] . The 
first factor is a little more subtle [33]: Considering that 
the transit time for a single particle through the slow 
site (in the absence of steric hindrance) is q^ 1 +£, q e ff is 
defined as the average rate to move just one step in this 
process: 



Qeff 



q£ 



l + ql 



The end result for the current is the solution to the alge- 
braic equation 



J = Qeff 



P+ (J) [1 - p_ (J)] 
£ - p- (J) £ 



Here, we give an explicit form (which displays the £ = 1 
limit well) 



Js 



where 



KL 



Q 



Q 



(1 - Q + VT - 2Q) £ 



(17) 



_ 2q£(l + q£) 
~~ (l + q + 2q£) : 



Note that, for any q < 1, this approach predicts that the 
current is less than the maximal value of \ 2 and so, the 
system is always in the H/L phase. 

The results of both mean-field predictions for J as 
a function of q arc shown in Fig. 14, along with two 
sets of data: J q (1) and J q (363). As expected, Jskl 
was purpose-built for two infinite TASEP's connected by 
a slow site and provides a better fit to the data with 
the blockage deep in the bulk (k = 363 in N = 1000). 
On the other hand, it is understandable that, e.g., for 
k = 1, a R must be very close to q. Thus, we may expect 
that the system will have maximal current for q > 1/ 
This behavior is confirmed by the data, as illustrated in 
the figure for £ = 12. Since Jnmf (q) has the property 
that it saturates at x 2 for q > x, it provides a better 
fit for J q (1). Of course, we recognize that, as mean- 
field theories, neither (14) nor (16) are the first step in a 
systematic expansion, so that they may better be thought 
of as "semi-phenomenological" . 
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FIG. 14: (Color online) Comparisons of the current, J, as a 
function of q. The legend labels the two sets of simulation 
data (slow site at k = 1 and 363) and predictions from two 
mean-field approximations. 



More importantly, there is a more serious, inherent 
limitation in this level of mean-field theory. It cannot 
account for the full fc-dependence in J q (k), since it deals 
only with infinite systems. One possibility to incorporate 
finite-size effects (of the left sublattice) in this kind of 
theory is to exploit the MG expression [7] 



Jmg 



PiPi+i 



P\+i 



for sites near the boundaries. Using it in a recur- 
sion relation for finding both the density profile and 
the current, we will arrive at a /c-dependent expression, 
Jmg = 1, fc), which replaces Jl in the first expres- 
sion of Eqn. (12). However, the high density fixed point 
of this recursion relation is unstable and very careful nu- 
merical analysis will be necessary. Preliminary work on 
this approach is promising and will be reported elsewhere 
[49]. 

We end this subsection by noting that the effects of 
a single defect in TASEP have also been investigated in 
[32]. Unlike our focus here - the dependence of J on 
the location of the slow site, they are concerned with a 
"multicritical system," i.e., a = f3 = 1/2 for the £ = 1 
case. Putting the defect at the center of the lattice, they 
explored density profiles in detail, finding power law tails 
on both sides of the defect with ^-dependent exponents. 
By contrast, our choice of a = = 1 places us far from 
the multicritical point. We have no reason to expect 
similar power laws. 



B. Two slow sites 

The most general TASEP with just two slow sites can 
be quite involved, since the parameter space is four- 
dimensional: {qi, <72, k%, k 2 }. To carry out a manageable 
investigation, we let both sites be deep in the bulk, so 
that only the distance between them, d = (k 2 — k\), plays 
a significant role. Further, as pointed out above, the cen- 
tral section resembles an ordinary TASEP with a and (3 
controlled by q\ and 52, respectively. Therefore, it is the 
smaller (slower) of the two rates which limits that cur- 
rent, which in turn dictates the current through the whole 
system. Thus, we will focus only on the q\ = q 2 = q case. 
Our parameter space will then resemble the single slow 
site case. 

Following the single defect case, the simplest levels of 
mean-field theory treat our system as three subsections 
with obvious labels: L, C, and R. From our discussion, 
only two (out of the many logical possibilities) combina- 
tions of phases, H/S/L and M/M/M, are expected to be 
viable. In addition to Eqn. (12), we have 



Jc = 



«c(l - etc) 
1 + a c l 



(18) 



Since the defect rates are identical, we fully expect that, 
for such a mean-field theory, (3c — etc- Now, matching 
the currents of the subsections, we immediately arrive at 
Jl = Jc = Jr and so, (3 L = a c = Pc = ct R . From 
here, the "naive" mean- field approach for H/S/L pro- 
ceeds identically to the above. The argument relies on 
the presence of a shock in the central section, so that 
there is a low (high) density region near site k\ + 1 (k 2 ) 
and we can impose the same discontinuity in the den- 
sities across both defects, i.e., p + (J) before and p- (J) 
after. Thus, we again arrive at Jnmf (q), given explicitly 
in Eqn. (15). The same argument can be applied to the 
next level of a mean-field approximation, which predicts 
Jskl (q), as in Eqn. (17). The major difference between 
the two approaches, as in the single slow site case, is the 
absence of the M/M/M phase in the latter. Meanwhile, 
their limitations are similar: The d dependence in J q (d) 
cannot be accommodated without serious modifications. 

Nevertheless, the spirit of these approximations can be 
exploited to provide J q (1) in the q — > limit. Since the 
central section consists of just one site, there can be no 
shock. Instead, the remnant of the shock is reflected in 
the average density there. It is more convenient to re- 
gard the system as two infinite TASEP 's, with nontrivial 
matching across a "doubly - slow site." Of course, we 
cannot expect to find any of the fascinating profile de- 
tails (e.g., inset of Fig. 10), but we should be able to 
obtain the "coarser" information, such as the currents. 
The goal is to understand the behavior of /S. 2 {q — > 0) 
(cf. Fig. 12) in, say, the first two non-vanishing orders 
in q. 

Now, for q <C 1, we are naturally in the H/L phase and 
the crudest approximation should suffice for the lowest 
order in the current. So, we let the bulk densities be 
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at their extremes (i.e., one and zero) and simply consider 
the time it takes for a particle to move through the block- 
age, from the moment its predecessor is "released." The 
current is just the inverse of this quantity, i.e., 



-+t-2 


-l 

_^ q 


A 


~* 2 



l-l(t-2) + ... 



(19) 



where we have included O (g 2 ) terms for computing the 
next order. But, at this order, we should also take into 
account that, occasionally the density before/after the 
blockage deviates from unity/zero by virtue of the right 
hand side of Eqn. (11) being non-zero. Thus, these den- 
sities are 

p L -> 1- J =l-q/2+ ... 
PR -» J£^q£/2 + ... 

and further suppress the current at the next-to-lowest 
order through the factor 



Mi-p«)^i-|(* + i) 

Combining these factors, we arrive at 



J^o (d = 1) - £ 



(20) 



l-9l<-2 



If we use exactly the same arguments for the q — > limit 
current in the case of one slow site, we find, instead of 
(19), 



1 



q[l-q{£-\) + ..] , (21) 



and, instead of (20), 

Pi(l-p«)-l-g(* + l) + ... • 

Finally since J q (d — > oo) is the same as the single- 
blockage current, we write 



J q ^o {d^oo)^q[l-2qe+...] 



(22) 



so that 



A 2 (g^0) 



It is remarkable how well this crude approximation 
agrees with the data in Fig. 12. There is no doubt 
that all curves extrapolate to the ^-independent value 
of 1/2 at q = 0. As for the slope at the origin, we can 
obtain a good estimate from the lowest q data points, 
using [A 2 (q = 0.02) - 0.5] /0.02. The values obtained 
from simulations for I = 1,2,4,6, and 12 are 0.92, 
1.55,2.53,3.44, and 6.06, respectively. 

We are aware that the expansion (22) differs from the 
small q limit of Jnmf- Unfortunately it is difficult to 



implement the same scheme for Jnmf here, 
must start from the exact pair of equations: 



since we 



J = l{fk (1 - n k +t)) = q(rk+i (1 - n k +e+i)) • (23) 

Various attempts at approximating p\ or p\ +1 led to 
poorer results. Alternatively we could exploit the ar- 
gument in SKL [33] and consider the average time to 
traverse both slow sites, 2/q + (1 — 2). This gives us a 
new effective q: 



Qeff 



ql 



2 + q(£-2) 



which can be inserted into Eqn. (16). The result is 
A 2 (q -» 0) -> \ + f (I + 2) + the O (q) term of which 
differs from the data by about a factor of 2. Clearly 
mean-field approaches are far from ideal for finding quan- 
titative predictions of J q (d). On the other hand, either 
Jnmf (q) or Jskl (?) provide tolerable results when the 
blockages are from from each other or the boundaries. 
Such variations in the quality of mean-field theories point 
to the importance of correlations. Considerable efforts 
appear to be necessary for a comprehensive, yet relatively 
simple, theory. 

Let us end this section with another method which 
could possibly improve the theoretical predictions [51], 
especially for the first few values of k. The idea is to 
find exact results, by solving the full master equation, for 
TASEPs (with extended objects) on very small lattices 
and then to match these to an infinite system (the R 
sublattice). One expectation is that, as in the I = 1 case, 
the finite-size current is larger than its counterpart for the 
infinite TASEP at the same (a, /3). This approach may 
eventually provide the essential argument to understand 
the increase in J q (k) as k becomes smaller. Similarly, 
this idea can be applied to the case with two slow sites 
when d is 0(1). Work is in progress to explore these 
avenues. 



SUMMARY AND CONCLUSIONS 



In this study we consider an inhomogeneous TASEP 
with open boundaries and populated with particles of fi- 
nite extent, I. The hopping rates arc uniform (set at 
unity) except for one or two sites ( "defect bonds" ) , where 
the rates, q, are different (faster or slower). We are inter- 
ested in the effects of these local defects on the density 
profiles and the currents through the system. Simulations 
with various t < 12 show that fast sites have no effect 
on the current, but induce discontinuities in the density 
profiles. In contrast, slow sites generate a significant re- 
duction of the current as well as nontrivial structures in 
the profiles, e.g., long tails behind the blockage, with pe- 
riod I. These findings are entirely consistent with similar 
studies in the past, most of which were restricted to £ = 1 
[31, 32]. If the inhomogencities are deep in the bulk and 
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far from each other, the current depends only on q and 
can be understood through simple mean-field considera- 
tions. Through the current-density relationship (for an 
infinite homogeneous TASEP), the overall densities in 
each of the defect-free sections can be also predicted, so 
that the various "phases" of these subsystems can be un- 
derstood. The distinguishing feature in our study is how 
the location of the defects affects the behavior of the sys- 
tem. For the case of one slow site, the current is slightly 
but measurably enhanced when the defect approaches 
the boundary. On the other hand, a drastic reduction of 
the current is observed when two slow sites are brought 
closer to each other. It is tempting to interpret these ef- 
fects as "interactions" between the defects and to seek a 
formulation that can describe them quantitatively. 

At present, neither the enhancement nor the suppres- 
sion of the currents can be understood in terms of simple 
mean-field theories. The essential limitation is that they 
are based on matching homogeneous TASEPs of infinite 
length. More sophisticated versions, relying on recursion 
relations for the particle density at each site, may be ex- 
ploited to deal with the finite subsections and provide 
some promise for a better understanding of such effects. 
For very small subsections, such as k, d < 5, it may be 
possible to find exact solutions (even for £ > 1) that 
can be used to match the mean-field descriptions for the 
macroscopic subsection(s). Work is in progress to inves- 
tigate these approaches systematically. 

Beyond one or two blockages, we should study systems 
with multiple slow sites, as our eventual goal is to under- 
stand the properties of fully inhomogeneous TASEPs. In 
particular, though restricted to just one or two slow sites, 
our findings - that the relative locations of blockages are 
important - will have implications for translation. For ex- 
ample, they are directly applicable to "designer genes" , 
which consist of many repeats of the same codon, except 
at one or two locations. Using the abundance of associ- 
ated aa-tRNAs as a control for the elongation rate across 



any particular codon, we can test our results directly on 
such genes. Thus, it will be interesting to see the physi- 
cal manifestation of, e.g., enhancement and suppression 
of production rates of such artificial proteins, depend- 
ing on the placement of the slow codon(s). Similarly, 
reproducing the intriguing ribosome density profiles will 
be revealing. More important than "designer genes," we 
should consider the implications for real genes. Our re- 
sults should provide, at the least, some simple qualitative 
insights. We can obviously maximize the production rate 
of a particular protein associated with a certain real gene 
by systematically replacing all slow codons with synony- 
mous, faster ones. However, for most genes this opera- 
tion will require a large number of replacements. Instead, 
with our findings in mind, we can achieve considerable 
increases in the production rates by making only a few 
substitutions, namely, by replacing the slowest codons, 
or a cluster of nearby slow codons. The ratio of cur- 
rent enhancement to the number of codon replacements 
may be used to quantify how "optimal" a certain set of 
substitutions is. This idea can be applied to finding op- 
timal means to suppress prodcution rates as well. Simu- 
lation work with TASEPs associated with real genes is in 
progress and we hope to demonstrate that these concepts 
are viable. 
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